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Abstract 


In recent years, powerful design tools for linear time-invariant multivariable control systems 
have been developed based on direct parameter optimization. In this report, an algorithm 
for reliable optimal control synthesis using parameter optimization is presented. Specifically, 
a robust numerical algorithm is developed for the evaluation of the // 2 -like cost functional 
and its gradients with respect to the controller design parameters. The method is specifi- 
cally designed to handle defective degenerate systems and is based on the well-known Pade 
series approximation of the matrix exponential. Numerical test problems in control synthe- 
sis for simple mechanical systems and for a flexible structure with densely packed modes 
illustrate positively the reliability of this method when compared to a method based on 
diagonalization. 

Several types of cost functions have been considered: a cost function for robust control 
consisting of a linear combination of quadratic objectives for deterministic and random 
disturbances, and one representing an upper bound on the quadratic objective for worst- 
case initial conditions. 

Finally, a framework for multivariable control synthesis has been developed combin- 
ing the concept of closed-loop transfer recovery with numerical parameter optimization. 
The procedure enables designers to synthesize not only observer-based controllers but also 
controllers of arbitrary order and structure. Numerical design solutions rely heavily on 
the robust algorithm due to the high order of the synthesis model and the presence of 
near-overlapping modes. The design approach is successfully applied to the design of a 
high-bandwidth control system for a rotorcraft. 
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Chapter 1 

Introduction 


1.1 Optimal Control Synthesis and Parameter Optimiza- 
tion 

Traditional design methods in linear optimal control for continuous-time systems have been 
extensively treated [1]. Development of these control systems is usually based on the char- 
acterization of the control problem under the setting of optimization of the two-norm of a 
set of controlled output responses to random disturbance inputs or initial conditions. Ad- 
ditional consideration of design robustness is taken by formulating the problem to include 
//“-norm bound constraints for a class of additive and multiplicative uncertainties applied 
at the plant inputs and/or outputs. Solutions are obtained for both the state- and output- 
feedback design problems and involve, in the majority of cases, solving an appropriate set 
of algebraic Riccati equations [4, 5]. Theoretical studies of these approaches have been a 
major concern of researchers in the control field and a major breakthrough has been made 
in recent work by Stoorvogel [13, 14]. 

An alternate and less often mentioned design option for robust multivariable control is 
based on direct numerical optimization of a H 2 performance objective with a given pre- 
specified controller structure. Early work in this area has been published by Levine and 
Athans [15], Anderson and Moore [16], and more recently, an extensive review of the subject 
was presented by Makila [17]. More recently, a new look into parameter optimization as it 
applies to multivariable control synthesis is provided by Ly [18] where he used a quadratic 
performance objective based on a finite-time horizon and the nonlinear optimization tech- 
nique of Ref. [21]. 

Synthesis of robust multivariable control systems using nonlinear constrained optimiza- 
tion has become increasingly important in recent years [40] [41] due to the design flexibility 
this approach offers. A numerical design package SANDY was developed to provide for 
the first time such an effective controller design tool. The plant and controller models are 
assumed to be linear and time invariant. The controller can possess a non-observer based 
structure and be of much lower order than the plant. To ensure good responses in criti- 
cal spectral regions, a frequency-shaped H 2 performance objective can be implemented by 
means of bandpass filters . Robustness can be enhanced by simultaneously optimizing over 
several plant model perturbations from a given nominal condition. Other control design 
specifications, including robustness, can be achieved via the concept of Closed-Loop Trans- 
fer Recovery [11]. This procedure, to be discussed further in Chapter 5, allows designers 
to achieve design performance and robustness starting from a satisfactory state-feedback 
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control law. 

Development of reliable numerical optimization for linear optimal control synthesis re- 
quires effective algorithms for the evaluation of the objective function and the associated 
gradients. One numerical limitation in the design optimization scheme developed in [18] 
arises when degenerate modes occur in the closed-loop system. The corresponding system 
dynamics matrix will contain defective eigenvalues which are most often visualized in Jordan 
blocks. Occurrence of these defective eigenvalues becomes more frequent when a multiloop 
control law is being synthesized based on frequency-shaped design objectives or distur- 
bances with degenerate power spectra (e.g. transverse Dryden turbulence spectra in wind 
turbulence). Degenerate systems can also occur in the design procedure of Closed-Loop 
Transfer Recovery (CLTR). Here, closed-loop dynamics from an output-feedback design in 
the CLTR procedure tend to overlap those attained under state feedback. With inexact 
arithmetic and large system matrices, these conditions lead to near degenerate systems and 
the appearance of Jordan blocks. 


1.2 Key Research Motivation 

• Experience with the optimization of low-order output feedback controllers in II 2 
inspires a search for a systematic way to include robustness in the design without 
tedious adjustment of design weights. One approach would be to use a different cost 
function, such as one based on worst-case initial conditions. Experience with this 
worst-case cost function indicates that while it would offer some help with robustness 
in particular, one would hope for a more comprehensive design methodology to obtain 
a more complete set of desired properties. Such is the promise of closed-loop transfer 
recovery (CLTR), which is based on recovering the properties of a state- feedback 
design with an output-feedback controller. 

• Numerical optimization based on gradient information promises fast convergence, but 
it requires a precise calculation of the cost function and its gradient with respect to 
the controller design parameters. Existing numerical algorithms for these calculations 
lose accuracy as degeneracy in a system mode is approached and several of the cor- 
responding eigenvalues become more defective. This loss of accuracy occurs in many 
design situations related to control of flexible structures with closely packed modes, or 
use of model matching and/or frequency-shaped objectives. Experience with closed- 
loop transfer recovery involving systems of reasonable size almost always generates 
defective eigenvalues. 

• The terminology “more defective” is not common usage. However, this terminology 
does properly express the quantitative nature of the loss of accuracy that prevails in 
existing eigenvalue-eigenvector decomposition routines. One can define a tolerance, 
based on the condition of the eigenvector matrix, where inaccuracies in the established 
gradient calculation will noticeably affect the optimization process. 

• In this research, several algorithms tolerant of defective systems have been considered. 
The most successful one involved converting the gradient computation to an exponen- 
tial calculation and using the well-known Pade series for the matrix exponential. The 
Pade series itself is applied to a scaled matrix, and the careful term-by-term evalua- 
tion combined with a judicious scaling process have enabled this approach to handle a 
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wide range of system matrices. In addition to the theoretical assurances, this formula 
was demonstrated to work in many selected test cases. 

• The proposed robust algorithm will provide an enabling technology for closed-loop 
transfer recovery within the framework of numerical optimization. System matri- 
ces under the CLTR design procedure tend to contain Jordan blocks owing to their 
size and the presence of overlapping poles. Within the CLTR framework, the best 
possible performance of an output-feedback controller can be predicted by exam- 
ining the system theoretic properties with the Special Coordinate Basis (SCB). In 
addition, a low-order output-feedback controller can be employed to recover many 
state-feedback properties such as robustness without the tedious effort of manipu- 
lating penalty weights and the many ensuing design iterations. A demonstration of 
CLTR on a rotorcraft problem containing a reasonably large system model is pre- 
sented. A systematic design procedure produces encouraging results over two flight 
conditions. The numerical optimization part of this procedure would not have been 
possible without the robust gradient formulation. 

• A more ideal approach for gradient calculation aims at combining robustness of the 
new gradient algorithm with the speed, low memory usage, and scale insensitivity of 
the original diagonalization method [18]. With this comes the need of an alternate 
means of decomposing a matrix into its eigenvalues and eigenvectors — with the priority 
of keeping the condition of its eigenvector matrix above a predetermined level for 
computational reliability. Such an algorithm is sketched in Appendices A and B. 


1.3 Summary of Contributions 

The basic contribution is a robust algorithm for computing the integrals 

X{t) = f \ As Be Cs ds , M{t)= f T z A{v ~ s) Be Cv De Es dsdv, 

Jo Jo Jo 

which are at the heart of the gradient calculations for several types of cost functions. There 
are two main developments: first, the ability to express each integral as a matrix exponential, 
allowing well-known techniques for computing a matrix exponential to come into play; 
second, the rearrangement of the matrix exponential computations into a faster, less memory 
intensive, and more robust form specific for each integral. 

Although the robust algorithm guarantees an answer for a wider range of input matrices, 
it uses more memory and is slower. The mode degeneracy is already parameterized in the 
condition of the eigenvector matrix. In running a wide variety of test cases, a limit on 
where this condition affects the optimization was determined. Using this limit, a means of 
switching between the original calculation and the robust one is made to promote optimum 
speed. 

Another contribution is in implementing an approximate minimax cost function: unlike 
// 2 -optimal control, this formulation minimizes a worst-case response to initial conditions. 
The actual worst case is approximated, and result of this approximation is in a form similar 
to that already used in // 2 -norm. 

Application of numerical optimization to CLTR completes a new development in CLTR 
design and analysis. A high-order and realistic design problem for the control of a UH-60 
rotorcraft poses a challenging limit on the number of optimization trials one may pursue due 
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to the CPU requirements. It is shown that, when properly set up, CLTR using numerical 
optimization would allow designers to obtain satisfactory design results quickly and in a 
routine manner for output-feedback controllers starting from a satisfactory state-feedback 
design. While the required controller order may be higher in the CLTR design approach 
than in a pure direct optimization of a H 2 - norm objective using an arbitrarily selected 
low-order controller structure, the need to re-design for robustness or other requirements 
is reduced. 
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Chapter 2 

Control Design Using Numerical 
Optimization 


2.1 Introduction 

Direct parameter optimization provides a versatile method for linear multivariable con- 
trol and has a broad range of applications. The design optimization is usually formulated 
within the context of H 2 optimal control. For completeness, a general formulation of the 
control design problem is given in this chapter. Different design cost functions are defined 
corresponding to a class of deterministic and stochastic control problems. Equivalent rela- 
tions between each design setup are established. Solutions of these design problems are a 
strong function of the disturbance input characteristics or plant initial conditions. These 
influences become particularly significant in the control synthesis of low-order controllers 
using direct parameter optimization. Such a design issue was addressed in [19] concerning 
the potential design sensitivity to plant initial conditions in optimal control synthesis. A 
worst-case design approach based on the largest singular value of a weighted covariance 
matrix was henceforth proposed by Bryson [19] to desensitize the optimal design. A sim- 
pler approach based on an upper bound to the worst-case cost is developed in this chapter 
and provides a convenient and numerically efficient way to address the worst-case design 
problem. A simple helicopter control design is used to illustrate the differences between 
a standard LQG design, an // 2 -optimal control using parameter optimization, and those 
obtained using different design techniques for worst-case initial conditions. 


2.2 Synthesis Model Description 

The following problem formulation is suited for the control synthesis of a robust low-order 
controller in linear time- invariant systems. The system P l (s) is to be controlled by a 
constant-gain controller C(s) as depicted in Figure 2.1 where z*(s) is the controlled output 
vector, y\{s) the measurement output vector, w l (s) the disturbance input vector and u'(s) 
the control input vector. For a consistent notation, the superscript i is used throughout to 
denote the system model at the i th plant condition. Although there may be a set of plant 
models corresponding to a series of design conditions, there is only one controller, C(s), in 
the design optimization. Hence, this single controller will provide stability and performance 
over all the defined plant conditions. The controller is modelled as a linear time-invariant 
system of a given pre-specified order. Its formulation can accomodate both a feedforward 
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Figure 2.1: Closed-Loop System with a Feedback/Feedforward Controller 

and a feedback controller structure. As stated before, plant parameter design robustness 
requirements in the context of our problem formulation are defined under the conditions 
that the control system C'(.s) stabilizes the class of plants P*(s), 1 < i < N p . 

State equations describing the system model P l (s) at the I th plant condition (Figure 2.1) 
are given in equations (2.1)-(2.3) below. 

• State Equations : 

/ X^t) = *V(<) + GV(f) + rV(f) 

UHO) = o [ZA) 

where x l {t) is a nxl plant state vector, u l (t) an mxl control vector, tu*(£) an m'x 1 disturbance- 
input vector, F 1 an nxn state matrix, G 1 an nxm control distribution matrix and P an nxm! 
input-disturbance distribution matrix. 

Notice that, in the above description, we assume that all the system states are initially 
acquiescent. This assumption is made without loss of generality since one can always 
establish impulsive inputs tc*(£) together with an appropriate influence matrix P to generate 
any given state initial conditions. 

• Measurement Equations : 

yl(t) = H\x\t) + D^u\t) + D^it) (2.2) 

where y l s (t) is a pxl measurement vector,//] a pxn state-to-measurement distribution matrix, 
D\ u a pxm cont rol- to-measurement distribution matrix and D\ w a pxm' input-disturbance- 
to-measurement distribution matrix. 

• Criterion Equations: 

z l {t) = + D l cuU l (t) + D i cw w i (t) (2.3) 

where z l {t) is a p'xl criterion vector, H l c a p'xn state-to-criterion distribution matrix, D 1 ^ a 
p'xm control- to-criterion distribution matrix and D 1 ^ a p'xm' input-disturbance-to-criterion 
distribution matrix. 


6 





For generality, the disturbances w l (t ) are modeled as outputs of a linear time-invariant 
system excited by either impulse inputs or white noise signals. In this manner, one can 
shape the disturbance signals to have any deterministic response (e.g. filtered step functions, 
sinusoidal functions, exponentially decayed or growing sinusoidal functions, etc...) or, to 
model stochastic inputs with any given power spectral density functions. At the i ih plant 
condition, the disturbance model is given by equations (2.4)-(2.5) below. 

• Disturbance State Equations : 

f xl(t) = Fi,xl(t) + T^rfit) , 4 v 

14 ,( 0 ) =0 M 

where x x w (t) is a n'xl disturbance state vector, rj l (t) a m'xl vector of either impulses (i.e., 
rf(t) = rf 0 6{t) with E[rfo] = 0 and E[r} l 0 r] l J'\ = W^), or white-noise processes 77(f) with zero 
mean and covariance E[rf 0 {t)ri x J {t)\ = W*6(t — r). The covariance matrix W l Q is an m'xm 1 
diagonal positive semi-definite matrix, F% an n'xn' state matrix of the disturbance model 
and r w an n'xm' input-distribution matrix. 

• Disturbance Output Equations : 

4(f) = Hi,x l w (t) + D i w rf i {t) (2.5) 

where w l (t) is a m'xl disturbance output vector, an m'xn' disturbance output matrix 
and D l w an m'xm ' direct feedthrough distribution matrix. 

State model of the controller C(s) (Figure 2.1) is that of a linear time-invariant system 
described by equations (2.6)-(2.7) below. 

• Controller State Equations : 

f x c (t) = A c x c (t ) + B c yl(t) . . 

Uc(o) = o 

where x c (t ) is a rxl controller state vector, A c a rxr state matrix of the controller and B c 
a rxp measurement-input distribution matrix. 


• Control Equations : 

u'(f) = C c x c (t) + D c y\(t) (2.7) 

where u l (t) is an mxl feedback control vector, C c an mxr control-output distribution matrix 
and D c an mxp direct feedthrough matrix. 

For control-law synthesis, any elements of the controller state matrices can be chosen as 
design parameters while the remaining elements can be left fixed at pre-assigned values. In 
addition, if needed, linear and nonlinear equality or inequality constraints can be established 
among the selected design parameters in order to ensure additional design structure. For 
convenience in the derivation of the performance index and its gradients with respect to the 
controller design parameters, we define a matrix C 0 that assembles all the controller state 
matrices (A c , B c , C c , D c ) into one compact form as follows, 


Co 


D c C c 
_B C A c 


1 (m+r)x(p+r) 


( 2 . 8 ) 


Thus, the single matrix C 0 will completely define the controller state model. Obviously, 
for the case of a static output-feedback design (i.e., the controller order r = 0), we simply 
have C 0 — D c . In the following sections, we will develop a control design problem based on 
the minimization of a performance objective using the controller C(s) defined in equations 
(2.6)-(2.7). 
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2.3 Formulation of the Closed-Loop System 

In the general case, the problem formulation assumes that there is no implicit-loop paths 
within the feedback control system. Namely, the control input u’(f) or the measurement 
output y\ (t) must not have any direct link to itself. This translates into the conditions that 
one of the products D c D l su or D l su D c must be zero. These conditions are not restrictive 
since, in practice, the presence of either actuation or sensor dynamics would automatically 
result in a system that satisfies the above assumptions. Moreover, for well-posed systems 
one can always reformulate the problem into an equivalent problem involving a modified 
set of measurement outputs j/*(f)(= y\{t) — D\ u u(t)). Thus, 

V l s(t) = HWit) + D i sw w i (t) (2.9) 

Note that the re-defined measurement outputs y\ (t) do not have a term involving the control 
variable u l (t) (i.e., D\ u = 0). However, for the discussions that follow, we make use of the 
results in [18] where we assume that either D l su D c = 0 or D C D ^ = 0. 


2.3.1 Case D c D l su = 0 


Combining the states of the plant, controller, and disturbances into an augmented system 
with the following states 

' At) 


At) = 


At) 

A(t) 


Dynamics for the overall closed-loop system are 

At) = F*At) + rV(t), 

where 


( 2 . 10 ) 


( 2 . 11 ) 


F n = 


F i + (PDcHl G i C c (P + G i D c D i sw )H i w 

B C {1 + D\ U D C )H\ A c + B c D\ u C c B c (\ + D\ u D c )Di w H' w 
OOP, 


r /l = 


(p + ^DcDDd* 
B c (l + Dl u D)Di w Dl 
P. 


w 
i i 

sw^w 


-I (n+r-|-n / )x(n+r-fn / ) 

( 2 . 12 ) 

(2.13) 


(n-f-r+n'JXm' 

+ Di u C c + 


px(n+r+n') 


Ds = 


(I + Dl u D c )Di w Dl 


pxm r 


II c = 


Hi + D^DcHl Dl u C c (D^DcD^ + D^Hi 


P / x(n-f-r-+-n / ) 


and 


(2.14) 

(2.15) 

(2.16) 

i , „ (2.17) 

With the above definitions, equations (2.2), (2.3), and (2.7) for y\(t), z'(t) and it* (t) become 

yi(t) = H'ix'^t) + (f) (2.18) 

At) = H'ix'\t) + {D^D c D^ + D^)DiAt) (2.19) 

At) = C H x ll (i) + D c Dl w Olr)'(t) (2.20) 


C H = 


D c Hl C c D c D\ w Hl 


8 



2.3.2 Case D l su D c = 0 

The closed-loop system for the combined plant, controller, and disturbance dynamics is 
again defined by 

x' i (t)=F ,i x' i (t) + r' i r j i (t) > (2.21) 

where the system matrix F H and the distribution matrix P are now given by 

F* + G i D c H i s CP(l + DcDiJCc (r + ^DcDiJHi’ 

B C H' S A c + B c D\ u C c B c D\ w Hl (2.22) 

0 0 F* 

' (Tt + GtDcDjjDi, ' 

T n = B c D\ w D l w . (2.23) 

r l 

Equations (2.18), (2.19), and (2.20) for the sensor, control, and criterion outputs are defined 
as before, where the matrices H' s \ D'j, H£, and C ,i take on the following form 

H'i = [Hi Dl u C c Di w H' w ] (2.24) 

D's = [D^Dl] (2.25) 

He = \H C + D^DeHs D^l + D c D su )C c (D CU D C D SW + D^H *] , (2.26) 

and 

C* = [l D C H\ (I + D c D' su )C c D c Di w Hl] . (2.27) 

2.4 Design Cost Functions for H 2 Optimal Control 

To examine the entire class of // 2 -optimal control problems and to handle the problem 
of sensitivity to plant modeling uncertainties, we define the objective function given in 
equations (2.28) and (2.32). This formulation turns out to be versatile and well-posed for 
the setting of a nonlinear constrained optimization problem. However, depending on the 
types of disturbance model, that is whether the disturbance outputs w*(t) are responses 
to impulse or white-noise inputs, different definitions of the objective function are needed. 
Another cost function is defined in equation (2.38) to address a worst-case control design 
problem. 

For these cost functions to be well-defined (bounded), one needs to identify the presence 
of direct feedthrough terms from the disturbance inputs ^(t) to the criterion outputs z'(f). 
Specifically, for either impulsive inputs or white-noise disturbances, the objective functions 
defined in equations (2.28), (2.32), and (2.38) would become unbounded when the direct 
feedthrough term D 1 ^ is nonzero. When a direct feedthrough term exists (e.g., in command- 
following synthesis), the disturbance inputs must be implemented as band-limited signals 
(i.e., they are generated from outputs of some roll-off shaping filters). Basically, this reduces 
the direct term to zero and places the direct influence of the disturbances within the 
criterion output distribution matrix H' l c in equation (2.16). 

A similar problem would occur for the case where there is a direct feedthrough term 
between the disturbances w(t) and the measurements y s (t). The condition for a bounded 
performance objective requires that the products Q i {D i cu D c D i sw + D‘J and R i D c D \ w be 
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zero for the cases of either impulsive or white-noise disturbances. Again, such conditions 
will always hold for band-limited sensor noises. 

These cost functions are given using the terminology defined in the previous section for 
the closed-loop systems. 


2.4.1 Random Impulsive Disturbances 

Consider the following cost function 

J\ (t/, Co) = W v r ^ 2 °‘ it E[z iT (t)Q i z i (t) + u iT (i)RV(t)]dt (2.28) 

2 , =1 Jo 

The expectation operator E[~] is over the ensemble of the random variables in the pa- 
rameterized impulse inputs rf{t) = rf- 0 6{t). Control design problems formulated with the 
above performance index J(tf,C 0 ) are often classified under the category of deterministic 
control. Under this category are, for example, the familiar control problems of command 
tracking control, disturbance rejection of unwanted but known external input signals, im- 
plicit and explicit model-following designs, // 2 -control to initial conditions and //^-control 
to sinusoidal inputs. 

This form for the cost function can be further expanded by embedding the closed-loop 
system responses 

x H (t) = (2-29) 

directly into the cost function C Q ) as 


Ji(t/,C 0 ) = 


- jr W* [ S E \rfj r * T e (F ' i+ 0 ,il)Tt [H ,i ?Q i H ,i c + C ,lT R l C H ]e (F '' +a ' !)Ti dlY ,l rf 0 

2 i=i Jo 

\ E W v E {€ Y nT S l (tf, Co) T'% 

1=1 


1=1 

AT, 


Co)T*E [rtrtj]} 

1 *=i 

1 n p 

2 E W P tr { r ' iT ^/. Co)r H w l 0 } 


(2.30) 


where 


S i (t f> C 0 )=[ f e (F '‘ +a</)Tt [R'c T Q i I/'c + (2.31) 

Jo 

and Wl = E [r&rtf 


2.4.2 Random White-Noise Disturbances 

The objective fimction defined in equation (2.28) is no longer valid since it is unbounded 
due to the presence of disturbances with white-noise spectra. For this reason, it is not 
proper to mix, in a given single plant condition, design considerations of initial conditions 
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and of disturbance rejection to white-noise disturbances with the same objective function. 
Alternate objective functions for the disturbance rejection problem must be defined (Refer 
to Appendix A of Ref. [18]). There are basically two ways to define a quadratic objective 
function for white-noise disturbances. The first formulation involves the state covariance 
responses evaluated at the terminal time tf. Namely, 


j *p 


Mtf,Co) = -^W^z^QV^) +u iT (t/)RV(f/)]. (2.32) 


i= 1 


The second formulation is a time- average of the state covariance responses, 


1 Np r rh 

Mtf,C 0 ) = —J2K E «' / * iT WQV(0 + « iT (f)RV(f) 

2 t J L/o 


dt 


(2.33) 


The expectation operator E a i[— ] is over the ensemble of the random processes defined in the 
input variables t f(t) for a closed-loop system destabilized by a factor a 1 . The destabilization 
effectively adds a value a 1 to the diagonal elements of the closed-loop system matrix. Given 
one of the above performance indices, one can address the entire class of // 2 -norm based 
control design problems. For example, we can solve for the linear quadratic regulator 
design (LQR), the linear quadratic gaussian (LQG) design, loop transfer recovery (LTR), 
closed-loop transfer recovery (CLTR), or model reduction based on the minimization of the 
// 2 -norm of the error. 

These two cost functions can be shown to be related to each other and also to the cost 
J\(tf,C 0 ) defined in Section 2.4.1. Embedding the closed-loop state responses, 

x H (t) = I*' e (^+«‘/)Ct- T )r 'y ( T ) dr, (2.34) 

Jo 


into C 0 ), we have 


N 

J 2 (tf,Co) = tr { ( H c lT Q lH c + C ,lT R l C ,i ) e (f’ , *+a*/)(t / -T) r /» 


E 


With E 


i~ 1 


N v 


(2.35) 


V(r)r? iT (s)] r ,l V F ' i+<W)T(t /- s) drds] 

7 ? t ( T ) r ? lT ( s )j = Wtf(r — s), equation (2.35) becomes 

N p 

Mt f ,D c ) = - Y^Wp tr { + C^FCC'^ 

J* f e (F>'+c l 'I)(L f -T) T n w i r nT e (F’'+ a i I) T (t f -T) j 
i W P tr { (Hc T Q IH c + C' iT R i C' i ^j e (F ,i +<* i nt T >i W ij : tiT e {F' i + a i J) T t dt 

2 i=i t Jq 


N 

= tr { r/iT f* e (F '‘ +aV)Tt ( HfQ i H * -I- C' iT R i C' i ') e (F ' <+a</)t dtT H W^\ 

t=i ^ ® ' 
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Or simplifying, 


Mtf, Dc) = \jr w; tr {r v&itf, D c )T*Wi] . (2.36) 

1= 1 

The above objective function is identical to the cost function D c ) in equation (2.30). 

As before, the objective function Mtf,C 0 ) can be simplified using the closed-loop re- 
sponses given in equation (2.34) as follows, 

N p 

Mtf , Co) = — Y, K tr { (h* t <?H* + C' iT R l C ,x ) *(*,)} (2.37) 


where 


*(«/) = Jo' LI e (F ' i+Qi/)(t - T) r 'E [^(t)^)] drdsdt 

Again, with E jto(r)u; T (s)] = W q 6{t — s ), 


Mtf, Co) = tr { ( h c T q 1 h ' c 1 + c ht r 1 c h ) 

Lo Jo e (F ' i+Q</)(t_r) r rt H^r rtT e (F ' <+ai/)7 ’( t - T > drdtj 

N 

= ~ J2 Wptr |r /tT £ £ e {F' i +a i 1 ) T {t-T) ( H *Tqi h H + C>iT R i C ,i j 

e (F'*+aV)(t-r) drdtT 'i W i^ 

— j-J q Mt,c 0 )dt < Mtf, Co) 


since the objective function Mt,C 0 ) is a monotonic increasing function of time t. Hence, 
Mtf ) Co) < Mtf, C 0 ) for any finite terminal time t /, and they become equal in the steady- 
state case where f/ — ► oo. 

We have shown that a design objective for random white-noise disturbances as given in 
the cost functional J2(tf,C 0 ) or J-s(t /, C 0 ) is equivalent to one involving initial conditions 
or impulsive disturbances. 


2.4.3 Worst-Case Impulsive Disturbances 

A worst-case design objective to impulsive disturbances, i.e. rf{t) = rf 0 8{t) is defined as 
follows, 


, JVp f tf [z ir (t)QV(f) + 

Mtf, Co) = - £ w; max -e 

z t=l % Vo Vo 

Substituting the closed-loop system responses given in equation (2.29), we have 


(2.38) 


Mtf,C 0 ) = -J2w; max 

A . - 71* 


i=l 


v l 7vl 
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The control problem for worst-case impulsive disturbances is to minimize the maximum 
eigenvalue of C 0 )r'\ or equivalently its maximum singular value since the matrix 

r nT S l (tf,C 0 )V l is symmetric and positive semi-definite. It is well-known that when singu- 
lar values of a matrix are repeated, they are not differentiable with respect to elements of 
the matrix ([23], p.288). With the relation a 2 (A ) < tr{A T A} for any arbitrary matrix A, 
we can re-define the worst-case functional by its upper bound, 

N p r 2 

Mtf,c 0 )<Mt f ,c 0 ) = Y,w l p tr{(r ,ir 5 i (i / ,C' 0 )r' i ) 

i=l ^ 

= j(s i (* / ,c 0 )r' i r' lT 

where S l (tf,C 0 ) is defined in equation (2.31). Derivatives of this upper bound are easy to 
compute. Thus, we prefer this approximation based on the Frobenius norm to the exact 
worst-case design objective. While it usually best to scale the disturbances through the T ri 
matrix, in practice, it is often convenient to allow an additional weighting matrix W a in the 
matrix r'T /lT as Y h WlY hT . Equation (2.39) becomes 



Jo 


( tf,Co ) =±w;tr{ (s*{ tf , c 0 )r' i w*r /iT ) 2 | 


(2.39) 


This corresponds to a bound on the maximum singular value problem 

JMI/.CW = | E K 

* i=1 nl Vo Vo 


The problem formulation for worst-case disturbances can also be extended to address 
cases where the disturbances have a given deterministic time behaviour. In this situation, 
the control problem is then to find an output-feedback controller for these specific types of 
disturbances (e.g., steps, ramps, sinusoidal inputs of known frequencies, and others). 


2.4.4 Worst-Case for White-Noise Disturbances 


The usual convention for white-noise disturbances is that the components to the disturbance 
vector are independent, and since E[rf(T)rf T (s)] = W’<5(r - s), this implies that W l a 
is full rank. However, for the worst-case design we optimize over W* that are rank 1, thus 
Wl — rforfj corresponds to a disturbance of the form rf(t) = rj l 0 P l (t), where p'(t) is a scalar 
white noise source of unit spectral density. Within this context we can write the two forms 
of the worst-case cost functional for white noise (analogous to the two forms of the H 2 cost 
for white noise) . The first is related to the state covariance response evaluated at the time 
tf 


1 N v 

J? c (tf, C 0 ) = - ^ Wp max 

z VV‘ 


i= 1 


tr{Wi} 


(2.40) 


The second is an average of the worst-case state covariance response in the time interval 

[o M 


j wc 


J t=l ° 


'[zV'WQiz*® + ?d T (0R I u l (t)] d 

tr{Wi} 


(2.41) 
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By embedding the dynamics, both cost functionals can be written in the form 


1* . tr{T HT S i (tf i C 0 )T H W*} 

Mt/, C 0 ) = -Y,K max — h 


i~ 1 


tr{Wi} 


The cost driven by such a noise corresponds to the maximum singular value of T' lT S l (t /, C 0 )T n 
and would be an upper bound to the regular H 2 cost defined for white noise. 

As before, we use the Frobenius norm of T llT S l (tf, C 0 )T' 1 to bound its maximum singular 
value. The result is the same as given in equation (2.39). 

Control-law synthesis using the above worst-case objectives and properties of their corre- 
sponding solutions are not well understood and need further investigation. This is a subject 
left for future research. However, it seems to have a potential of providing controller de- 
signs that are insensitive to plant disturbances in applications such as model reduction and 
closed-loop transfer recovery. 


2.4.5 Design Features of Direct Optimization 

Note that the performance indices given in equations (2.28), (2.32), and (2.38) are evaluated 
to a finite-time horizon tf. The use of a finite time plays a significant role in the imple- 
mentation of a reliable design algorithm for the optimum steady-state solution. It should 
be recognized that the objective function is well-defined regardless of whether the feedback 
control-law is stabilizing or not. Furthermore, a class of problems associated with command 
tracking of neutrally stable or unstable target responses (e.g. step and ramp commands, 
sinusoidal trajectories) is only tractable under the setting of a finite-time objective func- 
tion, but not in the confine of a steady-state objective function where tj — » oo. In practice, 
whenever possible, steady-state results are usually achieved when the terminal time tj is 
equal to five or six times the slowest time constant in the closed-loop system responses. 

Besides the concept of design based on a finite terminal time tj, other unique and 
important features have also been incorporated into the design objective function of equa- 
tions (2.28), (2.32), and (2.38). First of all, this objective function is not the usual quadratic 
cost function traditionally defined in linear optimal control. It is instead a weighted average 
of quadratic performance indices evaluated over a set of design conditions (1 < i < N p ). 
Different weights are assigned to each plant condition through the scalar variable W1 where 
W p > 0. Of course, design for only one plant condition (with N p = 1) reduces to the usual 
quadratic cost function evaluated at the nominal design condition. The time-weighted fac- 
tor e 2aH allows us to impose directly a stability constraint on the closed-loop eigenvalues in 
the // 2 -optimal control problems. Namely, when a steady-state design has been achieved 
and the optimum objective function is bounded, then the closed-loop system eigenvalues for 
the controllable modes will have real parts less than —a 1 . Finally, the weighting matrices Q l 
and R l are symmetric and positive semi-definite matrices. Note that our solution approach 
to the minimization of the objective function J(tj, C a ) is based on nonlinear optimization; 
hence, it does not require the control weighting matrix R 1 to be positive definite. In fact, 
in some design problems such as command tracking and model reduction, a proper objec- 
tive function contains only penalties on the tracking or model-matching errors and does 
not include penalties on the control variables (i.e., with R l = 0). The above observation 
further indicates the applicability of the design procedure in addressing the significant class 
of solvable singular optimal control problems [33]. 
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2.5 Derivation of Cost Function Gradients 


In this section, explicit gradient expressions for the cost functions defined in Section 2.4 with 
respect to the controller matrix C 0 are derived. Results are given for both cases D C D 3U = 0 
and D SU D C = 0. To begin, we consider the following useful differentiation rules of the trace 
function, 

-^-tr {QC o n} = {UQ) t and ~tr {QCjn} = KQ 

In addition, differentiation of the exponential function with respect to its argument can be 
obtained as follows. For a linear time-invariant system, 

4>{t) = A<t>(t ) 

m = i 


Then, 


d<f> ... . d<p . . dA 

dC 0 ® ~ A dCo ^ + dCo^ 


Solving for d(f>/dC 0> we obtain 

d<p 




dC< 

With the following definitions (as suggested in [18]), 


F l = 

X O 


F l 0 T l H l u 
0 0 0 
0 0 FL 


G'o = 


n = 


G { 0 

0 I 

0 0 


T l D l 


r 

1 w 


J (n+r+n')x(n+r+n') 


(n-j-r+n 7 ) x 


(n+r+n'jxm' 


H' = 


H l 0 D l H l 
0 I 0 

H\ = [Hi 0 D^Hi 


D l 0 = 


D\ = 


D\ 


D l D l 


J (p-fr)x(n+r+n') 
. p'x(rH-r+n') 

(p+r)xm' 


D'su o 

0 0 

= [dL 0 


(p+r)x(m+r) 

p'x(m-j-r) 


(2.42) 

(2.43) 

(2.44) 

(2.45) 

(2.46) 

(2.47) 

(2.48) 

(2.49) 


15 


To = 


n = 


T 1 - t 1 °]mx(m+r) 

0 0 
0 I 
0 0 


(n+r+n') X (m+r) 


0 0 0 
0 I 0 


(2.50) 

(2.51) 

(2.52) 


I (p+r')x (n+r+n') 

Gradients of the respective cost functions can then be expressed in a convenient form. 

2.5.1 Gradients of Ji(tf,C 0 ) (Case D SU D C — 0) 

Define 

F n = F 1 0 + (G i 0 + T 2 C 0 D\)C 0 H l 0 (2.53) 

r* = ri + (Gi + T 2 C 0 D\)C 0 Di (2.54) 

H'i = Hl + DiCoHi (2.55) 

C' 1 = T.CJIl (2.56) 


then 


dJ x 

dC a 


Nn 


C tf,Co ) = E K ( D 2 T QH'c + T^RC a )X\t f ,C 0 )H^ 


i— ] 


+(G i 0 + T 2 C 0 D\) T [M l (t f , C 0 )H ‘ T + C 0 )T n W l 0 D * T ] 

+T 2 r [M'(t f ,C 0 )Hi T + &(t Jl C 0 )T*W i 0 D?KE\C 0 f 


where 


X%,C 0 ) = f f e ^ i+ai ^ r /i Wj r ,iT e (F ' , + QfiI)Tt dt 
Jo 

S*( tf , Co) = f tf e (F ' ,+QM ) T t (H* T QH£ + C HT RC ,l )e( F ' i+a 'V t dt 
Jo 

Co) = f* f e (F '* + “‘ 1 )T(t_<T) (//^ T gH' i + C^RC H )e^ FH+a,iTit r H W*r 


e (F ' i+ail)Ta da dt 

2.5.2 Gradients of J\{tj,C 0 ) (Case D C D SU = 0) 
Define 


(2.57) 

(2.58) 

(2.59) 

(2.60) 


F ' 1 

= Fi + G i 0 C o (H t o + D\C 0 T 3 ) 

(2.61) 

p/i 

= r i + GiCoDi 

(2.62) 

H* 

= H\ + DiCo(H l 0 + D\C 0 T 3 ) 

(2.63) 

C'i 

= TxCoiHi + D\C 0 T 3 ) 

(2.64) 


then 


dJ 


Nr, 


dC 0 

■nT a ai 


htf, Co) = £>’ \{D?QH« + T?RC n )X'(t f , CoXH* + D\C 0 T^ 


1=1 


+ G l jM\t f) C 0 )[W 0 + D[C 0 T 3 } t + (G i 0 C 0 D\) r M l (t f ,C 0 )Tj + G l JS(t f ,C 0 )r n W 0 D? 

(2.65) 
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2.5.3 Gradients of C 0 ) 

The upper bound J a based on the Frobenius norm is a continuous and differentiable function 
of the controller design parameters. Computational effort in the gradients of Ja-(t/,C 0 ) is 
essentially the same as that for the cost function Ji(f/,C 0 ). Recall from equation (2.39) 
that 

Mtf, Co) = 2L K tr | (tfitf, Co)P W’r iT ) | (2.66) 

Next, we derive the gradients of C 0 ) for the cases D SU D C — 0 and D C D SU = 0. 

2.5.4 Case D SU D C = 0 

Using the definitions of F'\ P 1 , H n c , and C' 1 given in equations (2.53)-(2.56), we have 
dJ - ^ C A = £ w; (D?Q l H n c + Tj R'C' 1 ) Xtvsitf, C 0 )H'J 

aL, ° i=l 

+(C* + T 2 C 0 D\) T M'w s {t } ,C 0 )H? + T? C 0 )(D\C 0 H' 0 ) T 

+(G* + T 2 CoD\) T S{t } )T i W i 0 T iT S i {t }> C 0 )V ,i W 0 D^ 

+T f S\t { , C 0 )T i W l X T S i (t f , Co)r H W 0 {D\C 0 Di) T (2.67) 

where 

Xws(tf,c 0 ) = [ t A {F '' +Q ^ t r l w^r lT s i (t fy c 0 )r i wx T e^ F ' i+a ' r > Tt dt (2.68) 

JQ 

Mws(tf, Co) = £' J* e (^'*+°h) 7 ’( t-a) ^ H 'iTQi H >i + c' lT R l C'^ e (F ' <+ “‘ I)£ 

rw^sritf, c 0 )r i w' 0 r iT e( F ' i+ail)Ta dadt (2.69) 


2.5.5 Case D C D SU = 0 


Using the definitions of F'\ P*, H ll c , and C' { given in equations (2.61)-(2.64), we have 


dJojtf^Co) 

dC 0 


N P 

E W P { D ? Q iff c + n R l C H } X$y S (t fi C 0 )(IP 0 + D\C 0 T z ) t 
1=1 

+L»fcJ { 0*2 Q'H* + Tj WC' 1 } X' ws {t f , C 0 )Tj 

+Gl r 'M\y S (t f , C 0 )[Hi + D\C 0 T 3 ) t + (G' 0 C 0 D\) T M' ws (l f , C 0 )Tj 

+G iF s i (t J ,c 0 )r i w' 0 r iT s i (t f , c a )[ r* + (2.70) 


2.6 Special Design Problems 

For further study, we simplify the problem described in the previous section to a single 
nominal plant condition (i.e., N p = 1). The plant state model is given by equations (2.1)- 
(2.3) and has the following simplified form 


±{t) 

= Fx(t) + Gu(t ) 

(Plant dynamics) 

x(0) 

= x 0 

(Initial conditions) 

Vs(t) 

= H s x(t) 

(Measurement variables) 

z(t) 

= H c x(t ) 

(Criterion variables) 
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where we omit the superscript i (i.e., i = 1) for clarity. Note that simplification has been 
made to the measurement and criterion equations where direct influence from the control 
and disturbance inputs has been omitted (i.e, D su = D sw = = 0). For a 

deterministic control problem related to plant initial conditions, the disturbance inputs 
w(t) are not considered in the controller synthesis. Furthermore, it is assumed that the pair 
(F, G) is controllable, (F, H s ) is observable, and ( F,H C ) is detectable. 

In the conventional LQR problem, one seeks a control-law u(t) that minimizes the 
following cost functional 

J(t f , u(t )) = i J*' [ z T {t)Qz{t ) + u T {t)Ru(t)] dt (2.71) 

where the weighting matrices Q and R are both symmetric and positive definite. Solution to 
this dynamic optimization problem would normally proceed by adjoining the state dynamics 
x(t ) = Fx(t) + Gu(t ) as constraints to the cost functional J(tj) using a set of Lagrange 
multipliers X(t). From the stationary condition 6J(x,u, X,t/) = 0, we obtain the familiar 
Hamilton-Jacobi equations for a two-point boundary value problem. The resulting optimal 
control-law is given by u(t) = — R~ l G T X(t). By defining X(t) = S(t)x(t), and letting 
tf —* oo for the steady-state problem, one obtains a stabilizing static state- feedback control- 
law u(t ) = D c x(t ), where the gain matrix D c is given by D c = -R~ 1 G T S ss , and where S ss 
satisfies the algebraic Riccati equation 

F T S ss + S ss F - SssGR- 1 G t S ss + HjQH c = 0. (2.72) 


We now examine an alternative development to the above LQR problem in the per- 
spective of direct optimization using a fixed controller design structure. The problem is 
formulated as detailed in Section 2.4 and simplified to the case of a static output- feedback 
design (i.e., with C 0 = D c or u(t) — D c y s {t)). Clearly, the state-feedback case of LQR 
design is simply a special case where y s {t ) = Ix(t ), i.e., H s = I. Closed-loop system state 
responses to initial conditions under a static output-feedback law are given by 

x(t) = e^ GDcH ^x 0 = e p,t x 0 (2.73) 


The cost function J(tf,u(t)) in equation (2.71) can be rewritten as 

J(tj , D c ) = xlS{tj, D c )xo\ 


= ±t r {S(t f ,D c )E[x 0 xl}} (2.74) 

with S(tf , D c ) is given in equation (2.31), and the initial condition vector x 0 can be treated 
as a set of random variables with zero mean E[x 0 ] = 0 and covariance E[x 0 xJ} = W 0 > 0. 

The control design problem reduces then to the minimization of J(tf, D c ) with respect 
to the static output-feedback gain matrix D c . That is, 


dJ ( tf , D c ) 
dD c 


l-JL tr {S(t f> D c )W 0 } 


(2.75) 


Using Kleinman’s lemma [32] and after some manipulation, we obtain 


dJitf i D c ) 
dDc 


G t f 1 ' S(t f - a, D c )e( F,i + aiI >W 0 e( F ' i+a 'V T " daHj 
Jo 

+RD C H S f tf e (F i + a 'i)T Woe (F’'+o' i n T T dT H T ( 2 . 76 ) 

Jo 
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One could also obtain the above result using the general gradient expressions given in 
equations (2.57) or (2.65) when we specialize to the case of a single plant condition and 
C 0 = D c . 

Solution to the necessary condition for optimality 


dJ ( tf , D c ) 
dD c 


(2.77) 


for the optimal gain matrix D c is difficult when the terminal time tf is finite. However, for 
a stable closed-loop system matrix F' and as tf is increased to beyond five or six times the 
time constant of the slowest closed-loop modes, then 


S(t f - a, D c )e F '*W 0 e F ' T(J da - S(t f , D c ) j j* f e F '°W 0 e F ' T ° da j (2.78) 

Under this steady-state condition, we can re-assemble equation (2.76) as 

(rD c H s + G T S(t f , D c )) e {F ' i+aiI ^W 0 e( F ' i+a '^ T ° da) Hj ~ 0 (2.79) 

or, the term {hJ D jR + S(t f , D C )G ) is in the null space of the matrix 


H, f' 

Jo 


For the special case where H s = I (i.e., the state-feedback case) and W 0 > 0, we obtain the 
familiar result of 

D c ~ —R~ 1 G T S(t f , D c ) (2.80) 

To show the correspondence between the matrix S(tj, D c ) and the Riccati equation in the 
state-feedback case, we substitute equation (2.80) into equation (2.31) to obtain 

S(t f , D c ) ~ e (F-CR l GTs(t { ,D c ))Tt ( h t QHc + S ( tf> Dc )CR-'G T S(t f , D c )) 


e {F—GR~ 1 G T S(tf,D c ))t (2 gl) 

Multiplying the above equation by ( F T - S(t f , D c )GR~ l G T ) on the left and adding it to 
the same equation multiplied by (F - GR~ 1 G T S{t f> D c )) on the right, we obtain 



S(t f , D c )GR~ l G T ] S(t f , D c ) + S(t f , D c ) [F - GR~ 1 G T S(t f , D c ) 


,(F-GR- l G T S(t f ,D c )) T t 


HjQH c + S(t f , D c )GR~ x G T S{t s , D c ))}e^~ GR ~ lGTs ^ D ^ 


tf 


or 


o 

(2.82) 


F T S(t f , D c ) + S(t f , D C )F - S(t f , D c )GR~ 1 G T S(t f , D c ) + H F QH C = 

0 {F-GR~ 1 G T S{t s ,D c )) T tf 


[HjQH c + S(t f , D c )GR~ 1 G T S{t f , D c )\ e (F-GR-'G T s(t f ,D c ))t f 

(2.83) 

The above equation reduces to the algebraic Riccati equation under a stabilizing state 
feedback design when tf —> oo, 


F 1 S ss + S ss F - S ss GR~ l G T S ss + HjQH c - 0, 


(2.84) 
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where 


(2.85) 


S ss = lim S(t f , D c ). 

tf~+ OO 

Output feedback designs formulated in this manner yield the minimum cost functional 
value over an expected distribution of initial conditions. Often, as in the conventional 
state-feedback design, this expected distribution is assumed independent and uniform (i.e., 
E[x 0 xJ] = 7). 

In the static output-feedback design case (where H s ^ 7), analytical closed-form solu- 
tions cannot be found in general for the design problem stated in equation (2.28). One must 
resort to numerical techniques for the solution of the optimal design gain matrix D c . Note 
that, in the above formulation, solution uniqueness as well as the optimal solution itself will 
depend on the selection of the initial condition covariance matrix W a . Recently, Chen, et 
al. [33] have provided precise conditions under which an optimal static state-feedback design 
is unique. Non-uniqueness of an optimal state-feedback design can be used to explore other 
design issues such as robustness ([37]- [38]). 

One can also apply this special design case to the approximate worst-case cost. The 
stationary conditions given by 

m=m f{( s ^'^ rHf » rT ) 2 }=o 

become, for a sufficiently large t/, 

0 ~ ( RD C H S + G T S(t f , D e )) $*(t f )Hj (2.86) 

where 

&(t f ) = [ tf e {F ' i+aiI)a rw 0 r T s(tf, D c )rw 0 rV F ' <+0 ’ / > 7V da 

Jo 

One possible solution is given by the relation D C H S ~ -R~ 1 B T S(tf, D c ). Note that the 
above equation is identical to the H 2 form except for the weighting matrix 

{So e (F '’ +Q ' /)< TlF o r r S(f / , h t 

This is significant in fight of the fact that a dynamic output feedback problem can also be 
reformulated in terms of static output feedback [2] for Unear time-invariant systems, thus 
one can show that a similar weighting factor is responsible for the difference in H 2 and 
worst-case solutions for the dynamic output feedback case as well. Clearly, this solution 
is identical to the solution of the LQR problem when H s = 7. Generally, solutions of the 
above stationary condition for the optimal static gain matrix D c can be only be done using 
numerical optimization. 

2.7 Design Example Using a Simplified Helicopter Model 

A comparison of designs achieved under the LQ, 77 2 , and worst -case performance index is 
done for the case of state feedback. The state feedback result also establishes a baseline for 
the output feedback designs. 

For output feedback, we compare designs from three different performance objectives: 
LQG, H 2 , and worst-case. Because the disturbance set-up in the synthesis model is so 
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important, this comparison is performed over 2 sets of synthesis models. The first set is 
related to loop transfer recovery (LTR). The second set involves introducing fictitious noises 
into each of the states of the plant model with different intensity levels. This is an attempt 
to approach state-feedback results by improving the stability of the closed-loop system. 

2.7.1 Helicopter Model 

We consider a 1-dimensional model of the OH6A helicopter in hover mode (see Ref. [19]) 
with a longitudinal cyclic control 6(t). The vehicle states are pitch rate q(t), pitch angle 9(t), 
ground velocity u(t), and position x(t). The system units are feet, seconds, and degrees. 
The state model is 

X u X q -g 0 u(t ) ’ X s 

M u Mq 0 0 q(t) Ms r/,', , 

0 1 0 0 9{t) + 0 ° w + 

1 0 0 0 J [ x(t) J [ 0 

Values of the constants in the model matrices are X u = —0.0257, X q = 0.013, M u = 1.26, 
M q = -1.765, X s = 0.086, and M 6 = -7.4080. 

The vehicle is disturbed by a horizontal gust u w (t). It is modeled as a white noise with 
a spectral density of 18 ft 2 /sec and entered in all designs as well as in the calculation the 
RMS gust responses. 

2.7.2 Full— State Feedback Design 

Initially we consider a full state-feedback controller to minimize the following quadratic cost 
objective 

J(tf,D c ) = jf 7 E [; x 2 (t ) + £ 2 (f)j dt (2.88) 

which penalizes the ground position x(t) and the longitudinal cyclic control 6(t) equally. 
The optimum design D c is found to be D c = [1.9890, -0.2560, -0.7589, 1.00]. This solution 
can be obtained both from the Riccati equation (2.84) and from direct optimization based 
on equation (2.76). In the gradient search procedure, the function nelder in MATLAB 
was used to minimize the following objective functions: the objective function J(tj, D c ) in 
equation (2.88), the maximum singular value J<, for the exact worst-case design defined in 
equation (2.38), and the Frobenius norm for the approximate worst-case design in equa- 
tion (2.39). In all numerical cases, the terminal time tf is set equal to 512 seconds, which 
meets the aforementioned criterion for being steady-state. 


Table 2.1 

Static State Feedback for the Helicopter Model 


Eigenvalues 

Single-Loop Margins 

Unit Gust 

Open-Loop 

Closed-Loop 

Loop Phase Gain 

Responses 

-1.8891 

-1.8461 

x : 60.01° 8.61dB 

x = 0.243 

0 

-1.1192 

9: 41.4° -5.85dB 

9 = 1.497 

0.0492 ± 0.4608i 

(C=-0.106) 

-0.4464 ± 0.9774i 

«=0.415) 

6: 60.2° -11.3dB 

6 = 0.976 


Design results are shown in Table 2.1. They provide the baseline for comparison with 
the output feedback designs that follow. Time responses to a step command in ground 



" u(t) ' 

m 
. »(*) . 
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position, x CTnd , are shown in Figure 2.2 along with the Bode magnitude plot of the transfer 
function x(s) / x C md(s) . 

Design solutions are identical in all three cases and are equal to those obtained from 
the Riccati equation for the LQ problem. Note again that the LQ design is equivalent to 
the minimization of an objective function evaluated to a uniformly distributed set of initial 
conditions on all four rotorcraft states, i.e., with E[x 0 x„\ — I. In the numerical optimiza- 
tion, the initial conditions enter into the evaluation of the different types of cost functions 
explicitly. In the minimization of tr{S(tf, D c )} for LQ and H 2 problems, a(S(tj , D c )) in 
worst-case, and tr{S 2 (t/, D c )} in the approximate worst case, the maximum singular value 
of S(tf, D c ) is found to be 5.58. 

2.7.3 Comparison of Output-Feedback Designs 

To make a valid comparison, the controller order must be the same for each of the LQG, 
H 2 , and worst case designs. Since the LQG design results in a 4th order controller with 
no direct feedthrough terms, all other designs based on direct optimization will involve a 
4th order controller with D c set to zero. 

The LQG formulation is set as to the form of the disturbance inputs it uses, though 
being an observer-based design, it actually marries two different disturbance models. The 
state- feedback gains were determined using a unit uniform uncorrelated distribution on the 
initial conditions of all the rotorcraft states. The full-order observer part of the LQG design 
assumes white noise for both process and sensors. The H 2 and worst-case problems will 
also be formulated with white noise. We define the sensor outputs for the rotorcraft as the 
rotorcraft position x and pitch angle 9. Sensor noise in the ground position measurement 
x(t) has a spectral density of 0.4 ft 2 sec, and for the pitch measurement 9(t) a sensor noise 
of 0.4 deg 2 sec. The helicopter is excited by a horizontal gust of spectral density 18 ft 2 /sec. 

An initial output-feedback design is performed using only the gust disturbance as pro- 
cess noise. Compared to the results in Table 2.1 for state feedback, the RMS responses 
to the gust input are higher in all three output-feedback designs (Table 2.2). Robustness 
evaluation is conducted at the control actuator input. The resulting gain and phase mar- 
gins are low. The maximum singular value of S in the LQG design is somewhat larger 
than that for the state feedback design, and those from the H 2 and worst-case designs 
are significantly higher still. The approximate worst -case design does indeed optimize to a 
smaller a-(r r ST) than the other designs. This trend is maintained for the different sets of 
disturbances. Note that “tuning” of the worst-case design to gust disturbances leads to a 
poor step response characteristic. (Figure 2.4) Responses to command inputs for the LQG 
design always matches corresponding responses of the LQ design. 

2.7.4 Actuator Disturbance Augmentation 

Here we augment the disturbance model to include a process noise entering into the control 
actuator by letting T = [r a aG] . The power spectral density of this process noise is fixed 
at 18 deg 2 sec. (Note that the gust disturbance power spectral density is 18 ft 2 /sec.) The 
factor a in the disturbance distribution matrix is used to scale the effects of this additional 
process noise relative to the gust disturbance, and has the sequence of values [0,0.1,1,10]. 
This process simply follows the traditional Loop Transfer Recovery (LTR) procedure in 
LQG designs, to enhance the robustness in the actuator loop [39]. Table 2.2 shows the 
expected improvement in robustness. Responses to the gust disturbance also improved 
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Table 2.2: Traditional Loop Transfer Recovery Using Different Design Methods 
LQG Designs 


Fictitious Actuator 
Noise Factor a 

0.0 

0.1 

1.0 

10.0 

Mean Square x: 

1.255 

1.165 

0.627 

0.431 

Responses 6\ 

6.617 

6.125 

3.034 

1.888 

to iLm <5 . 

1.274 

1.254 

1.069 

0.950 

x: P.M. (deg) 

31.66 

31.56 

30.83 

30.51 

G.M. (db) 

5.04 

5.02 

4.90 

4.88 

6: P.M. (deg) 

23.93 

24.93 

30.30 

30.09 

G.M. (db) 

10.63,-4.110 

10.95,-4.117 

15.63,-3.569 

-3.75 

6: P.M. (deg) 

29.84 

29.90 

40.62 

52.49 

G.M. (db) 

10.66,-7.46 

10.97,-7.638 

15.62,-9.265 

-10.55 

a(S): 

28.57 

28.81 

30.66 

31.56 

c r(r ,T ST'): 

3.882 

4.442 

53.48 

4111 


H 2 Numerical Designs 


Fictitious 
Actuator Noise 

0.0 

0.1 

1.0 

10.0 

Mean Square x: 

1.276 

1.182 

0.627 

0.431 

Responses 0: 

6.596 

6.102 

3.034 

1.889 

to u w 6 : 

1.269 

1.244 

1.069 

0.950 

x: P.M. (deg) 

32.35 

33.31 

30.83 

30.51 

G.M. (db) 

5.15 

5.15 

4.91 

4.88 

0: P.M. (deg) 

24.23 

24.75 

30.30 

30.09 

G.M. (db) 

10.45,-3.268 

10.75,-3.296 

15.61,-3.569 

-3.75 

6: P.M. (deg) 

27.97 

28.46 

40.58 

52.47 

G.M. (db) 

10.49,-7.86 

10.78,-7.940 

15.60,-9.260 

-10.55 

a(S): 

12138 

20527 

2263. 

725.3 

a(r' i £T): 

3.889 

4.447 

53.52 

4112 


Worst-Case Numerical Designs 


Fictitious 
Actuator Noise 

0.0 

0.1 

1.0 

10.0 

Mean Square x: 

0.948 

0.887 

0.529 

0.342 

Responses 9: 

6.887 

6.209 

2.798 

1.901 

to Uw 6 : 

1.398 

1.343 

1.036 

0.965 

x: RM. (deg) 

29.24 

29.78 

30.08 

30.72 

G.M. (db) 

4.03 

4.01 

4.10 

5.31 

6: P.M. (deg) 

24.75 

23.56 

29.40 

30.36 

G.M. (db) 

11.46,-2.828 

12.03,-2.785 

18.86,-2.981 

-3.57 

6: P.M. (deg) 

29.67 

27.53 

43.37 

55.00 

G.M. (db) 

11.49,-6.18 

12.04,-6.514 

18.87,-9.495 

-10.80 

a(S): 

9448. 

16894 

1840. 

122.1 

6 f(r ,T ST'): 

3.562 

4.072 

50.40 

3972 
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when we increase the actuator process noise. Although there is no general tendency for 
a(S) to converge toward the state-feedback result, this value did decrease somewhat with 
increasing a in the H 2 and worst-case designs. For different values of a, the LQG and H 2 
designs show significantly different results in a(S) and in step responses to £crmi- 


2.7.5 Fictitious State Noise Augmentation 

In this set of designs, we augment the disturbance model with independent fictitious noises 
entering into each of the plant states by defining T — [r g al]. We set the power spectral 
density of each fictitious noise at 18. Again the scaling factor a in the disturbance distribu- 
tion matrix takes on a sequence of values [0,0.1,1,10]. The ensuing stability augmentation 
improves the overall controller performance, with the parameter a(S), representative of 
the state-feedback performance, trending toward the state-feedback value of 5.58 as a is 
increased. 

In the LQG design, increasing the process noises in each state will cause the estimated 
states of the observer to converge more rapidly. As a result, robustness at the actuator 
loop worsens, while the maximum singular value of 5 approaches that of the state feedback 
design (Table 2.3). The single-loop robustness at the sensors improves, and the controller 
is less attuned to gust alleviation in the position x (as in the II 2 based controller designs). 
The step responses and the transfer functions x/Xcmd are equal to those of state feedback 
as expected. 

Results based on the numerical optimization of H 2 costs do not correspond to those of 
LQG as well as in the previous case, though overall they have similar gust responses and 
robustness margins. Increasing the relative emphasis of fictitious disturbances going into 
the plant states caused no convergence of H 2 and LQG design properties. 

The step responses in the II 2 designs (Figure 2.5) converge nearly to the state feedback 
response as a is increased, though with a slightly slower response and less overshoot. While 
a(S) became smaller, did not settle to the state feedback value. 

For the worst-case designs, increasing the fictitious plant state noises slightly improved 
the overall disturbance rejection, while the robustness for all loops improved more substan- 
tially. Of note is that the robustness at the actuator loop improved. The value of a(S) 
tends towards that of the state feedback case, while the step responses (Figure 2.6) also 
converge to nearly the state feedback response. However, in spite of increased system sta- 
bility, disturbing the initial conditions of the states did not result in as desirable overall 
controller characteristics as augmenting with actuator disturbances. More favorable results 
may be obtained with larger disturbance amplitudes. 


2.8 Conclusions 

Direct parameter optimization provides a powerful means to address a wide class of design 
criteria and controller structures. The control formulation presented in this chapter shows 
several common types that are well-suited for the framework of numerical optimization. 
Development of reliable numerical optimization techniques is tied closely to the availability 
of reliable computation of the cost functions and its gradients with respect to the controller 
design variables. In all cases, evaluation of the cost functions presented in this chapter 
and their gradients involves directly the computation of the following integrals of matrix 
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Table 2.3: Singular Value Type Disturbance Augmentation 

LQG Designs 


Fictitious State 
Noise Factor a 

0.0 

0.1 

1.0 

10.0 

Mean Square x: 

1.255 

1.410 

2.653 

2.674 

Responses 6: 

6.617 

6.769 

9.130 

8.913 

to 8: 

1.274 

1.215 

1.072 

1.017 

x: P.M. (deg) 

31.66 

30.96 

34.77 

35.76 

G.M. (db) 

5.04 

4.44 

5.83 

6.58 

0: P.M. (deg) 

23.93 

22.76 

25.60 

28.03 

G.M. (db) 

10.63,-4.110 

10.91,-2.806 

12.83,-3.686 

15.66,-4.06 

6: P.M. (deg) 

29.84 

27.87 

25.30 

25.60 

G.M. (db) 

10.66,-7.460 

10.97,-8.10 

12.84,-8.516 

14.41,-8.63 

a(sy. 

28.57 

16.36 

8.450 

7.553 


3.882 

5.647 

126.0 

11140 " 


H 2 Numerical Designs 


Fictitious State 
Noise Factor a 

0.0 

0.1 

1.0 

10.0 

Mean Square x: 

1.276 

1.488 

2.656 

2.053 

Responses 0: 

6.596 

7.082 

9.049 

7.810 

to u w 8: 

1.269 

1.187 

1.064 

1.040 

x: P.M. (deg) 

32.35 

32.56 

34.81 

34.84 

G.M. (db) 

5.15 

4.66 

5.83 

6.20 

0: P.M. (deg) 

24.23 

23.64 

25.82 

27.79 

G.M. (db) 

10.45,-3.268 

11.25,-2.839 

12.90,-3.698 

16.22, -3.7^ 

6: P.M. (deg) 

27.97 

28.83 

25.58 

26.65 

G.M. (db) 

10.49,-7.86 

11.33,-7.703 

12.40,-8.562 

16.02, -8.6C 

a(S): 

12137 

1038. 

4300. 

70.11 

d(r ,i 5r'): 

3.889 

5.690 

126.1 

11101 


Worst- Case Numerical Designs 


Fictitious State 
Noise Factor a 

0.0 

0.1 

1.0 

10.0 

Mean Square x: 

0.948 

0.777 

1.597 

1.432 

Responses 0: 

6.887 

5.714 

6.720 

6.090 

t O Wy; <5 . 

1.398 

1.403 

1.049 

0.967 

x: P.M. (deg) 

29.24 

29.20 

32.71 

33.33 

G.M. (db) 

4.03 

3.96 

5.96 

6.26 

0: P.M. (deg) 

24.75 

22.97 

28.69 

30.93 

G.M. (db) 

11.46,-2.828 

11.40,-2.623 

16.82,-3.637 

-3.79 

6: P.M. (deg) 

29.67 

27.79 

28.88 

29.86 

G.M. (db) 

11.49,-6.18 

11.43,-6.977 

16.70,-8.596 

-8.87 

<t(S): 

9448. 

6770. 

296.0 

9.523 

<7(r ,7 sr'): 

3.562 

5.221 

113.1 

10413 
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exponentials 


X(t)= [ l e Ar Be° T dr 

Jo 

M(t)= [ f e A{v ~ s) Be Cv De Es ds < 
Jo Jo 


Note that the matrix S(t) has the same form as X(t). Efficient algorithms have already been 
developed for the above integrals based on diagonalization of the system matrix. However, 
they are prone to inaccuracies and lead to convergence problems in numerical optimization 
when the system matrix contains defective degenerate modes. A reliable algorithm for 
evaluating these integrals has been developed and its details are given in Section 3.3. 
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Figure 2.2: State-Feedback Design for the Helicopter Model 
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Figure 2.5: H 2 Designs for Various Stability Weights 
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Chapter 3 

Evaluating X(t) and M(t) 


3.1 Introduction 

Established methods for evaluating X(t) and M{t) are based on diagonalization of the 
system matrices A, C and E in the exponential functions. Diagonalization is achieved 
using a similarity transformation derived from the eigenvalue-eigenvector decomposition of 
these matrices. It is further assumed that the similarity transformation constructed from 
the eigenvector matrix is nonsingular. Let V A , Vc and V E be the eigenvector matrices of 
the matrix A, C and E respectively, then 

A = V A \ A VX l , C = V c AcVc l , E = V E A E V^, (3.1) 

where A^, Ac and \ E are diagonal matrices. One can express the exponential function of 

At 

e ™ as 

e At = e V A \ A VX l t _ y Ae A A ty-\ ( 3 . 2 ) 

Application of this decomposition in the calculation of X(t) is shown below. 

X(t) = f l e Ar Be° T dr 
Jo 

= V A jjf e A ^ T 6e AcT dr| V^ 1 (3.3) 

where B = V A l BVc • Advantage of this approach is that the exponential function of a diag- 
onal matrix is also diagonal. In this case, time integration in X(t) can be performed directly 
by explicit integration of a product of scalar exponential functions. The resulting numerical 
algorithm is quite accurate and efficient, provided that the transformation matrices V A and 
Vc are not ill-conditioned. A similar procedure can also be applied to the evaluation of 
M(t). Complete discussion can be found in Appendices C and D of [18]. Breakdown of 
this algorithm will occur when the matrices A or C develop Jordan blocks. This situation 
happens frequently in the control synthesis of flexible structures with densely packed modes 
as demonstrated in the design example of Section 4.2. 

Clearly, in order to have a reliable design algorithm for optimal low-order output- 
feedback control synthesis [18], one must develop a robust numerical scheme to evaluate 
matrix integrals of the form shown in X(t) and M{t) in the case where the system has 
defective eigenvalues. 
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3.2 Alternative Approaches for Solving X(t ) and M(t) 

One simple approach is to evaluate 

X(t) = j\ Ar Be CT dr , M(t) = J* Be Cv De Es ds dv (3.4) 

directly using a numerical quadrature. Efficiency of numerical integration techniques is 
poor; especially when it requires small integration step size for satisfactory accuracy in 
the case of stiff system matrices A, C and E. Another possibility is to use some types 
of algebraic Lyapunov equations for the solution of X(t) and M(t). For example, it can 
be easily shown that the matrix X(t) can be obtained from the solution of the following 
Lyapunov equation, 

AX(t) + X(t)C= [e 4r £e CT ] ‘ (3.5) 

Solution of equation (3.5) exists if A* (A) + Xj{C) ^ 0. As A, (A) + Aj(C) tends toward zero, 
the solution accuracy will degrade on a continuous basis. Furthermore, there are many 
situations encountered in practice where this scheme will run into difficulty. For example, 
when A = C T and the system matrix has poles at the origin. Thus, for practical purposes 
X(t) cannot be solved from a scheme based on Lyapunov equations. 

The Lyapunov equation for M(t) is even more poorly behaved. By changing the order 
of integration, the double integral in M 

M = f r e A[v ~ s) Be Cv De Es dsdv 
Jo Jo 

can be re-written as 

M = J* e A{q - r) Be~ Cr dr j e ct De E{t ~ q) dq 

Because the matrices A , B , and C are time invariant, one can rewrite the convolution 
integral as 

J e Ar Be c< ^~^ dr = [I 0] exp j ^ j 

Applying this to the inner integral of M, we obtain 

M = lo {[ I0 1 GXP { 0 -C I | e Ct De E ( t - q) dq 

This form is used to derive the Lyapunov equation for M. 

With the above equation, we write 

ME = j(|[I0]exp| g [ I ||« C ‘ d « E,e * 

= | [i o] ex P {[ o -c } <* - ' «)} [ i } } eC ‘ DeE, \ 

+ i*{ |101 0 -c exp { 0 -c <i_,) } I }‘ aD ^ E <^ 
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+ 


0 - \ [I 0] exp I 
[A B] exp 


0 


A B 
0 -C 

A B 
0 -C 


}[:]} 


(t~q) 


e Ct Dl 


e ct De Eq E dq 


Using the relation exp ■ 


A B 
0 -C 


t) = 


[ e 4£ /o e A ^~ T ^Be—Cr dr 


-Ct 


, we have 


ME = - [{^V (t - r) Be- Cr dr}e Ct £>] + Be~ c ^~ q) e Ct De Eq dq 

+ Jo A {Jo 9 em ~ Q) ~ r)Be ~ Cr dr } e Ct De Eq dq 


= B J* e~ c{t - q) e ct De Eq dq - {j* e Mt ~ r) Be~ Cr rfrj e Ct D 


+a£ |[I 0] exp 1 1 ^ B 


c 


(t-q) 


:} 


e Ct De Eq E dq 


Finally, we can assemble the Lyapunov equation for Ai as 

ME - AM = B j l e~ c{t ~ q) e Ct De Eq dq - \^j\ A{t ~ r) Be~ Cr drj e Ct D 

In most applications, the matrix E is equal to A and hence renders the Lyapunov equation 
for M unsolvable. 


3.3 Matrix Exponential Approach 


Another possible approach is based on the direct use of the exponential of a matrix. It is 
well-known [23] that convolution integrals involving matrix exponentials, as represented in 
the matrices X(t) and M(t), can be derived from the matrix exponential of an augmented 
matrix. It can be shown that the matrix X (i) can be derived from the upper-right partition 
of the following matrix exponential, 


X{t) = e M [I 




(3.6) 


Thus, computation of X ( t ) now involves the computation of a matrix exponential. A reliable 
algorithm for computing the matrix exponential is given in Section 3.4. 

In a similar fashion, one can express the matrix M(t) in terms of a submatrix of a 
matrix exponential. To see this, we start from its definition 

M(t) = f 1 f e A(v ~ s) Be Cv De Es dsdv 
Jo Jo 

e Av BeCvdv^ De Es ds 

Jo e ~ AS [Jt eAVBeCV dv } DeES ds 



35 


Let’s perform a change of integration variable v = t- r. We have, 


M(t) = j\- As y^ S e A ^- r) Be c ^dr^De Es ds 

= ~ It eA[t ~ S) {Jo 3 e ~ ArBe ~° r dr } e ct De~ E ^ d(t - s ) e Et 

e~ Ar Be~ Cr dr | e ct De~ Eq dqe Et 

= £ ^£ e A(q ~ r) Be Ct/2 e~ Cr drj e ct/2 De E{l ~ q) dq (3.8) 

Notice that part of the integrand in equation (3.7) delimited by braces can be replaced by 
terms involving the exponential of an augmented matrix. This follows directly from results 
developed for the matrix X(t). With this substitution, we obtain 

M(t) = j[/0]expj q Be _£ gj J | e ct/2 De E{t ~ q) dq 

= J q [I0] exp | q Bt _£ (t - g)| ^ j e ct / 2 l£ e Eq dq 

' [ A Be Ct ! 2 0 1 1 [ 0 

= [/ 0 0] exp i 0 -C e ct ! 2 D 1 1 0 (3.9) 

0 0 El 

V L J L J 

In this section, we have shown that the matrices X{t) and M(t) can be formulated in terms 
of solutions of matrix exponentials. Their evaluation depends therefore strongly on the 
accuracy and reliability of numerical methods for computing matrix exponential. We will 
present one such algorithm in Section 3.4. However for computational expediency, special 
consideration must also be taken to ensure efficiency of the overall scheme when the upper 
limit t is large and one of the matrices A, C or D is unstable. Also one must economize 
memory requirements associated with high dimensionality of the augmented matrix when 
computing the matrix exponential. These considerations will be elaborated in Sections 3.6 
and 3.7 where we give precise algorithms for the computation of the matrices X(t) and 
M(t) respectively. 

3.4 Numerical Method for the Matrix Exponential 

Several numerical methods are available for the computation of the matrix exponential [22]. 
Among these, an approximation method based on Pade series is found to be satisfactory [23]. 
An important component in any numerical routine for matrix exponential is the scaling of 
the matrix argument prior to the series calculation. Due to the simple result that e At = 
(e At / 2 ) 2 , a scale factor in terms of powers of two (i.e. 2 m ) is often used. In this scheme, one 
can recover the actual value of the original matrix exponential by performing m squarings 
on the matrix exponential of the scaled matrix. The index m is determined based on the 
desired size for the scaled matrix. In our algorithm, scaling is applied to the original matrix 
until its oo-norm [| A)^ falls below 1/2. 

As mentioned above, the preferred series approximation in our computation of the matrix 
exponential is the Padd series. Let’s review some of the unique features associated with the 
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Pad6 series for the case of a scalar function E(z). On its most basic terms, it is a rational 
function of z of a preselected order that approximates the function F(z). For a given choice 
of the order of the numerator (say TV ) and of the denominator (say M), the Taylor series 
representation of this Pad4 series must match the power series representation of T(z) for 
the first (TV + M + 1) terms. Namely, 


m ~ PmW = 


£,=o 


(3.10) 


In fact, the most common form of the Pade series is known as the diagonal sequence where 
the numerator and the denominator have the same order (i.e. M = N). While it is known 
that the Pade series for the matrix exponential E(z) = e z converges only slightly faster than 
the Taylor series for a scalar argument, the improvement becomes significant for a matrix 
argument. In the matrix case, Pade series involves computation of a numerator matrix 
Af(At) and of a denominator matrix V(At). For a diagonal Pade series of order TV , we have 


A f(At) = 1 + 




(27V -1)! TV! 
(27V)! (TV — 1)! 

(2 AT - i)! TV! 


(27V — 2)! TV! 


(2N)\i\(N -i)\ 


(27V)! 2! (TV 
(Aty + 


+ 


2 )\ 

TV! 

(2 TV)! 


(Aty 


(At) 


N 


and 


V(At) 


(27V — 1)! TV! (27V — 2)! TV! 2 

(27V)! (TV - 1)! (27V)! 2! (TV- 2) ! 1 ^ ’ 


+ (-l) 1 


(2TV — i)! TV! 
(2TV)!i! (TV - i)! 


(Aty H 1- (- 




The matrix exponential is simply given by 

e At = T>~ 1 (At) M (At) 


(3.11) 


Invertibility of V(At) is ensured by proper scaling of the matrix argument At. 

Another important consideration in the Pade series is its length TV. Assuming that 
the matrix At has been scaled such that is less than 1/2, the parameter TV can be 

choosen according to [23] such that 


o3-2JV (TV!) 2 < 

(2 TV)! (2 TV + 1)! ~ 


(3.12) 


where e is a given desired tolerance for accuracy. 

With TV determined in the above manner, the nominal error in a Pade series approxi- 
mation can be thought of as the exact calculation of a matrix exponential for a “nearby” 
matrix (At + E) where E is the error matrix with ||E||oo < e||AT||oo- The relative error of 
the approximation is bounded by the following inequality, 


e {At+E) _ e At 


< c||At|| 00 e £ ll j4t ll° 0 


(3.13) 


Thus, reducing the oo-norm of the matrix At would indeed improve the numerical accuracy 
of the matrix exponential. It has also been shown that methods by series approximation 
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yield better accuracy if the matrix argument has been preconditioned [29]. Additional 
improvement may therefore be gained by first preconditioning the original matrix. Another 
immediate benefit of lowering the oo-norm of the matrix being exponentiated is that the 
actual scaling factor m needed would also be smaller; thereby resulting in a fewer number of 
matrix multiplications in the squaring procedure. As usual, preconditioning a matrix tends 
to bring singular values of that matrix closer together (i.e. lower the condition number), 
thus avoiding situation where the scaling factor is predominantly determined by a few large 
singular values, and causing significant loss of precision related to the set of small singular 
values. The most common method used in the precondition of a matrix is Osborne’s method 
[25], which minimizes the Frobenius norm of that matrix (and thus indirectly lowering its 
oo-norm). However, extensive tests conducted so far seem to indicate that preconditioning 
of a matrix does not yield significant reduction in the oo-norm and a smaller scaling factor. 
However, Osborne’s method is relatively light in terms of computational burden. The 
procedure of preconditioning a matrix is nonetheless recommended from the point of view 
of improved accuracy ([26], [29]). 

In the implementation of our design algorithm for optimal low-order controller synthesis 
[18], a value of e = 10 -8 has been selected, requiring therefore a 4-term Pade series (i.e., 
AT = 4) in the evaluation of the matrices A' i (f), C(t) and M i (t) of equations (2.58)-(2.60). 
Additional considerations in the implementation of the proposed method for computing 
X{t) and M{t) are given in Sections 3.6 and 3.7. 


3.5 Preconditioning With Scaling and Rotation 

We can liken the original cost function calculation method (that of an eigenvalue-eigenvector 
decomposition) as a method that does best when the eigenvalues are widely scattered (at 
least, distinct). The matrix series calculation excels when a matrix is homogenized. For the 
particular use in a controller optimization scheme, the matrices encountered usually have 
widely spaced eigenvalues, but degeneracies in these matrices may still exist. 

Scaling the argument matrix to a Pade series to render its oo-norm less than 1/2 ac- 
commodates a large magnitude eigenvalue (stable or unstable). This same scaling would 
lead to roundoff errors for a small magnitude eigenvalue in the same matrix. This potential 
problem renewed interest in finding some method that went beyond Osborne’s method in 
homogenizing a matrix. We attempt such a method by combining the row and column 
scaling of Osborne’s method with an inverse Jacobi rotation. 

Suppose that one chooses a pair of columns and their corresponding rows in a matrix. 
Can one find a single rotation and two scaling terms to reduce the Frobcn ins norm of the 
matrix? Consider 


A i+ i = T~ l AiTi 


cos 8 

sin B 

<T1 

<?2 

sin 6 

cos $ 

v 1 

C2 


An ' Aij 
Aji • ■ • Ajj 
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X 


(t i cos 9 <7 1 sin 9 


—(72 sin 9 u 2 cos 9 


we want to determine a rotation angle 9 and two scaling terms o\ and <72 such that 


min 

0 , ( 7 - 1 , 0- 2 


T-'ATi 


It would then be hoped that a succession of these transformations would lead A n toward a 
minimum Frobenius norm value. The Frobenius norm was chosen to see if an easy means 
for determining 9 , <7i, and <72 was even possible. Unfortunately, there is no easy means for 
finding these scalings or rotations. The amount of work needed is significant, and thus this 
approach was abandoned. It seemed that, among all the alternatives, Osborne’s method 
was still the best. 


3.6 Detailed Algorithm for Computing X(t) 


As seen in Section 3.3, the matrix X(t) can be evaluated in terms of a matrix exponential 
as shown in equation (3.6). Conceptually, it is a simple and straightforward procedure to 
compute the matrix exponential of any arbitrary matrix using the Pade series discussed 
in Section 3.4. However, it becomes a nontrivial task when we try to implement an effi- 
cient algorithm that examines carefully the issues related to accuracy, speed and memory 
requirements. To understand the basic difficulties, consider in detail the components of the 
matrix exponential used according to our problem defined in equations (2.58) and (2.59) 
for the matrices X l (tf) and £*(£/) 


exp 



' -A B' 


■ e -At e -Atx ( t ) ' 


1 

O 

o 

l 

0 e Ct 


(3.14) 


where A — C T = F n +a l I. Clearly, if the system matrix A is stable (i.e. all the eigenvalues 
of the matrix A have negative real parts), one could easily encounter numerical overflow 
when evaluating the term e~ At even though the matrix integrals X(t) and C(t) are perfectly 
well-behaved. The overflow problem occurs most likely in the final squaring process. To 
arrive at a feasible approach in the evaluation of X(t), one needs to examine in detail the 
steps taken in arriving at the matrix exponential of the original matrix starting from that 
of a scaled matrix (i.e. in the squaring process). 

Let’s assume that one has scaled the input matrix A by AAt where At is a reasonably 
small time interval given by A t = t/n — t/2 m . Thus, we need to first evaluate 


exp 


- A B 
0 C 



where At = t/n = t/2 m . 


For notation convenience, we define 


exp 


-A B ' 

Atj = 

’ e~ AAt e~ AAt / 0 At e AT Be CT dr ' 


' D 

E ' 

0 C 

0 e CAt 


0 

F 


(3.15) 
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Furthermore, let W = exp(AAt) = D l . Now we can write our result as follows, 


X(t) = W n [D n ~ l E + D n ~ 2 EF + D n ~ 3 EF 2 + • • • + EF n ~ 1 }, (3.16) 


or 


X{t) = W[E + WEF + W 2 EF 2 + • • • + W n ~ l EF n ~ x \. 


(3.17) 


The above results are produced by performing m squarings of 


and taking the 


D E 
0 F 

2.58)-(2.59)), the solution 


appropriate submatrix for X{t). In our application (cf. equations 
would therefore involve products of matrices of size 2(n + r + n'). Close examination of 
equation (3.17) leads to the following algorithm involving only product of matrices of size 
(n + r + n ') with the final result achievable in m steps, 


Step 1 : 

Pi 

- W, 

Q i 

= E, 

R i 

= F 

Step 2 : 

P2 

= Pi 

Q 2 

= Qi + PiQi-Ri, 

R 2 

= R\ 



= w 2 , 


= E + WEF, 


= F 2 

Step m : 

Pm 

’ ’ ' ) 

— p2 

~ * 771—1 > 

Qm 

__ ' ’ j 

~ Qm— 1 5 

Rm 

= Rl 



= W n , 


+Pi m—\Qm—\ Rm— 1 > 


= F n 


Finally, X(t) = WQ m . It should be noted that one can ’’absorb” this extra factor of W 
(_ e AAt^ j n t G m atrix Q\ without any change to the above algorithm (i.e. starting the 
above algorithm with Q\ = WE instead). This removes the need to retain the matrix W 
throughout the computation. 

Finally, one notes that the terms Pi or Ri for (i — l,m) may underflow and become a 
null matrix for some z ; in particular when the scaling factor is large (i.e. m large). When 
this situation happens, one can simply truncate the series calculation for X(t) up to the i th 
step in the above algorithm since all of the significant (and nonzero) terms have already 
been accumulated into the matrix Qi. 


3.7 Detailed Algorithm for Computing M(t) 


Here the numerical algorithm is a bit involved compared to the one given for X(t). This is 
largely due to the increased complexity of the argument of the matrix exponential. Following 
the procedure described in Section 3.6, let’s perform a scaling upon the input matrix A by 
A At such that computation of the matrix quantities M 0 , H,J,P,U and W = V~ l is well- 
behaved. These quantities are defined from the following matrix exponential, 


r 

' A 

Be ct/2 

0 

1 


‘ P 

He ct ! 2 

Mo 


0 

-C 

e ct / 2 o 

At 

► = 

0 

V 

e Ct/2j 

k 

0 

0 

E 

J 


0 

0 

u 


(3.18) 


Due to the possible numerical underflow in the matrix e Ct / 2 for large t, the matrices H and 
J are computed directly from the following definitions, 

/•At 

H= e Ar Be Cr dTe~ CAt (3.19) 

Jo 
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and 


/•At 

j = e - CAt e CT De ET dr (3.20) 

Jo 

However, the computation of M 0 in equation (3.18) can still underflow due to its explicit 
dependence on e ct / ’ 2 . For the calculation of the matrix M(t), ideally it can be obtained from 
m squarings of equation (3.18). If carried out in this manner, potential numerical overflow 
is eminent since, according to our equation for M l {tj) in (2.60), we have A T = C = E T = 
F* + a l I. Hence, if the matrix C is stable, then the matrix exponential e~ Ct = V n will 
become unbounded. To bypass this difficulty, as in the calculation for X(t), one needs to 
conduct the squaring algorithm explicitly. It can be shown that the matrix M(t) can be 
computed as 


M(t) = P n -'M 0 + P n ~ 2 M Q U + P n ~ 3 M 0 U 2 + • • • + PM 0 U n ~ 2 + M 0 U n - } 


+ HW 2 J + HW 3 JU + PHW 3 J + PHW*JU + P 2 HW a J + • • • + P n ~ 2 HW n J (3.21) 


This formulation no longer involves the matrix V. The above series for M can be decom- 
posed into two parts — one that contains the matrix M 0 and the other that does not. The 
terms involving M a can be thought of as 


[7 0 ] 


’ p 

Mo ' 

n 

" o ' 

0 

U 


7 


( 3 . 22 ) 


which can be performed by m squarings. The remaining terms involving 77, J, IF, P, and 
U are computationally intensive and are of the form 


n — 2 n— 2 

EE P i HW 2+i + 5 JW where 2 + z + j < n. 

i— 0 j— 0 


(3.23) 


Without the restriction that 2 + i + j < n, this would have been formed as the product of 
two easily computed series, 


(77 + PHW + P 2 HW 2 + • • -)W 2 (J + WJU + W 2 JU 2 + ■■■) (3.24) 


An efficient procedure for computing the final matrix M.{t) is to merge both the easily 
computed portion given in equation (3.22) and the more difficult series in equation (3.23) 
into a sequence of m steps: 
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Step 0 : 
M 0 


H 

J 

p 

u 

w 

Step 1 : 

Mi = PM 0 + M 0 U 


H\ = H 

Ji — J 

p 2 

u 2 

w 2 

+HW n J 


+PHW 

+WJU 




Step 2 : 

m 2 = P 2 M 1 +M 1 U 2 


Ih = H ! 

J 2 = J\ 

pi 

t/ 4 

w 4 

+H\W n ~ 2 J\ 


+P 2 HiW 2 

+W 2 J l U 2 




Step 3 : 

m 3 = p 4 m 2 + m 2 u 4 


II 

lSj 1 

J 3 = J 2 

p 8 

u 8 

w 8 

+H 2 W n ~ 6 J 2 


+p 4 h 2 w 4 

+W i J 2 U i 




Step j : 

Mj = F 2 ’ 1 Mj-i + Mj-i 

u 2 *- 1 

Hj = Hj-i 

II 

pV 

u v 

w 2 ’ 



+P*' 1 Hj^W 2 *' 1 

+w 2i 1 Jj_ it/ 2 ’" 




Step m : 

Mm = P n/2 M m -i + M m - 

-1 lW 2 

P 771 = U 771 1 

Jm ‘ J m — 1 

pn 

u n 

w n 

+ 

1 

to 

1 


+P n / 2 H m - X W n ! 2 

-\-W n t 2 J m -\U n / 2 





where the final matrix M{t) = M m + H m W 2 J m . Due to potential numerical underflow, 
the term W l ~ 2 is not accurately obtained from the product W l V 2 where V — W~ l . Indeed, 
one needs to recompute the term W n+2 ~ 2,3 at each step of the above algorithm. This could 
become the major drawback in our scheme even though we have used an efficient matrix 
exponentiation routine that computes W 1 requiring at most 2log2{i) matrix multiplies. 

Further simplification of the above algorithm can be achieved if we make use of the fact 
that we have A = E = C T (cf. equation (2.60)) and therefore U — P = W T . If in addition 
W n is zero (or effectively so), restriction on the indices i and j of 2 + i + j < n in equation 
(3.23) becomes inconsequential. Hence we can express 


n~ 2 n— 2 

EE P i HW 2+i+j JU j = (H + PHW + P 2 HW 2 + • • -)W 2 (J + WJU -I- W 2 JU 2 + ■ • •) 

1 = 0 j= 0 


resulting in a simpler algorithm that involves the following three series, 


' p 

Mo ' 

71 

' 0 ’ 

0 
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I 


(b) ( H + PHW + P 2 HW 2 + • • •) 

(c) (J + WJU + W 2 JU 2 + ■■■) 
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This algorithm can again be computed in m steps as follows, 


Step 0 : 
Mo 

H 

J 

Step 1 : 

M\ — P M 0 + M 0 U 

Hi = H + PHW 

J x = J + WJU 

Step 2 : 

M 2 = P 2 M l +M 1 U 2 

h 2 = h x + p 2 h 1 w 2 

J 2 = Ji + W 2 J X U 2 

Step 3 : 

m 3 = p*m 2 + m 2 u 4 

H 3 = H 2 + P^lhw* 

J 3 = J 2 + W i J 2 U i 


Step j : 


Mi = P 2 i ~ l M<-i 

3 

II 

1 

Jj — Jj—\ 

+ 

Vj. 

1 

to 

Vi, 

+P 2 * Hj-xW 2 *" 

+W 2 i ~' Jj-xU 2 ^ 1 

Step m : 

M m = P n/ 2 M m - 1 

Hm ~~ F m 1 

1m = Jm — 1 

+M m -xU n/2 

+P n / 2 H m ^xW n/2 

+W n l 2 J m _xU n/2 


If, for some index j ^ m, W l (and likewise U l and P l ) is zero or nearly so, then the 
calculation of M(t) is reduced to M{t) = HjW 2 Jj since Mj = 0 where i = 2 J . 


3.8 Basic Test Results 

A direct evaluation of X and Ad on a degenerate system matrix provides a clear under- 
standing of the severity of the errors in the diagonalization method. Consider 


F = 


T = 


H'jQlV c + C^RC' = 

where F is constructed under a similarity 
transformation T, 
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transformation from a core matrix A and a 



1.0 

1.0 

-1.0 

T = 

2.0 

-2.0 
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-3.0 

3.0 
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Namely, F = TAT 1 . In the X computation, the results were 


' 58.37 

58.97 

30.58 ' 


' 61.90 

66.04 

19.98 ' 

58.97 

57.83 

45.91 

5 ^ robust — 

66.04 

71.96 

24.70 

30.58 

45.90 

00.00 
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24.70 
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As one can see, there is really no correspondence between the two results. Furthermore, in 
comparing the results for M, we have 


Adding — 


-3.623 x 10 9 
3.081 x 10 17 
2.054 x 10 17 


4.997 x 10 9 
6.162 x 10 17 
4.108 x 10 17 


-7.495 x 10 9 
-9.244 x 10 17 , 

-6.162 x 10 17 


•M-robust 


2211. 

2519. 

1230. 

1860. 

2146. 

1153. 

1817. 2133. 

1475. 


With these large discrepancies, it is clear that design optimization based on the diagonal- 
ization method would fail when defective eigenvalues appear in the closed-loop system. 


3.9 A Two-Mass- Spring Design Problem 

In this model, the open-loop system has defective eigenvalues at the origin; thus success 
of the diagonalization-based approach depends on the removal of this degeneracy through 
feedback. Hence, selection of the initial controller guess plays a key role. We consider two 
setups to illustrate the problems encountered by the algorithm based on diagonalization. 



Figure 3.1: A Two-Mass-Spring Mass System 

Equations for the dynamic model are given below. 

wii il\ — k(y 2 — Vi) + u + w 
m 2 y '2 = k(yi - y 2 ) 
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(3.25) 


(3.26) 


wucic nti — ii t 2 — m — *2 — i. me pro oiem is to control tne displacement or the second 
mass by applying a force to the first mass as shown in Figure 3.1. At the start, it is simple 
to verify that the basic open-loop system has a pair of defective eigenvalues at the origin. 
We have chosen a second-order controllable canonical form for our controller model, 


A = 


0 1 
A 22 


B = 


[C n C 12 J; 


D = \D n 
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The control design problem is to minimize the // 2 -norm of the following closed-loop transfer 
function T y2W between the disturbance w and the displacement 3/2 of the second mass using 
the controller design parameters -4 2 i , A 22 , C u , C X2 and Du. From the initial design guesses, 

^21 = —2, A 22 = — 1) C\\ = 0, C \2 — 0.5 and Du = 0 

The original method based on diagonalization in SANDY [18] converges to the optimal 
design gains, 

A 21 = —0.8571, A 22 = — 0.9258, C\\ = 0, C12 = —0.4535 and D u = —0.2449 

and the optimum value of ||T^ 2UJ ||§ = 7.71838215122. A summary of the resulting closed-loop 
eigenvalues is given in Table 3.1. 

Table 3. 1 : Closed-Loop Eigenvalues of the Two-Mass-Spri ng System 


Eigenvalue 

Damping 

Freq (Hz) 

-0.2290 

± 

0.3397i 

0.559 

0.0652 

-0.1553 

± 

0.8480i 

0.180 

0.1372 

-0.0786 

± 

1.29501 

0.061 

0.2065 


One cannot initiate the search using a compensator design of zero gains (i.e., A21 = 
A 22 = C n = C12 = Du = 0) because, in this case, the closed-loop system would have two 
pairs of degenerate eigenvalues at the origin; one for the rigid-body mode and the other from 
the open-loop compensator poles. The early algorithm based on diagonalization recognizes 
that Jordan blocks are present and the numerical search is halted. Detection of a Jordan 
block within SANDY is done by looking at the condition of the eigenvector matrix. If this 
condition goes above 10 12 or its inverse is less than 10“ 12 , the design search is terminated 
with a failed status. Condition of the eigenvector matrix influences the accuracy of the 
matrices X and M. At some point, the loss in accuracy will affect the overall optimization 
process. To alleviate this problem, it was suggested that one simply re-start the design 
algorithm with any compensator design (stabilizing or not) that initially produces a non- 
degenerate closed-loop system. This approach does not provide a satisfactory solution in 
general. 

A more insidious possibility is exemplified by the following initial design guess, 

-421 = —2, A 22 - — 1, C11 — 0, C12 — 0 and Du = 0 

It turns out that this case of Jordan blocks is not detected from the condition of the 
eigenvector matrix. The diagonalization-based algorithm fluttered around the initial point 
and finally terminated with no improvement to the original design guess. The condition 
number of the eigenvector matrix associated with this problem was 2.302 x 10 6 , which would 
generally be considered acceptable. 

The robust form of the algorithm allows convergence of the controller gains in all the 
three initial guesses shown. The final result is the same in each case and equal to 

A 2 i = -0.8571, A 22 = -0.9258, C n = 0,C X2 = -0.4535 and D n = -0.2449. 

The result is the same as the one determined from a successful run using the diagonalization 
algorithm. A history of the inverse of the condition number of the eigenvector matrix for 
the initial design guess of 

A 21 = —2, -4 22 = — 1, Cn = 0,Ci2 = 0.5, and Dn = 0, 
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Figure 3.2: Behavior of Condition Number over the Entire Design Optimization 


^60 


is shown in Figure 3.2. One can see that condition of the eigenvector matrix varies signifi- 
cantly over the entire run. Conceivably, the inverse condition number may jump back to a 
low value in a typical design iteration. Not very visible in this plot is the value at the first 
iteration of 1.8 x 10 -8 for the chosen initial guesses. After the first iteration, the inverse 
condition number jumps to 0.166. 

The main difference between the robust (Pade-series-based) form and the diagonalized 
form is in the CPU time of the overall computation. Results are obtained for a VAX/VMS- 
Workstation DEC-3500 as follows: CPU time of 19.59 seconds with the algorithm based 
on diagonalization, and 97.36 seconds using the proposed robust algorithm. The increase 
in computational load is expected and constitutes the basic trade-off between reliability 
and speed. The proposed algorithm is more reliable and with this advantage does take 
somewhat longer in computational time. 

Although a simple fix to the diagonalized form is to start with a non-degenerate closed- 
loop system by using a different initial compensator design, it has been found that the 
algorithm could occasionally break down due to the presence of near degeneracies. Thus, 
for a reliable design method, the solution algorithm must treat degeneracies as a common 
occurrence. This situation is more evident in the optimal output-feedback control design for 
high-order structural models with closedly packed flexible modes, and in the design proce- 
dure of closed-loop transfer recovery. These cases include a situation where two nominally 
independent modes appear as a Jordan block due to roundoff error in a higher order system 
(usually 30 -order or higher). 

For efficiency, the numerical algorithm currently implemented combines the robust al- 
gorithm with the faster diagonal form through a “switch” based on a threshold of 2 x 1 0 ~ 5 
in the inverse condition of the eigenvector matrix. This value is found to be adequate in 
detecting defective degenerate cases and computational accuracy of the run is improved by 
the intermittent calls to the robust form. 
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3.10 Degeneracy During Optimization 


Presented here is a model reduction problem that develops a Jordan block during the 
optimization process [31]. Consider a model with closely spaced modes, 



-1 

0 

0 ' 


’ l.i ' 

F = 

0.0005 

-1.000001 

0 

, r = 

1.2 


0.0005 

0.0005 

-1.00001 


1.3 


D = \ 0 ] 


Parameter optimization is used to reduce this model to a second-order system by minimizing 
the H 2 - norm of the error. Starting with the following arbitrarily chosen guess 


Ac 


-0.9995 1 

2 -1 


B r = 


-1.414 
1.414 ’ 


C c = [-5.616 0] , D c = [0] 


The optimization converges to the following reduced-order model, 


’ -1.027 

0.012 

r> 

’ -1.440 ' 

1.436 

-1.625 

j £> — 

1.542 


C = [-5.627 -0.107], 


D=l 0 ] 


with a final cost function of 4.76 x 10~ 7 . Diagnostics associated with this run indicate that 
while the gains were “optimized”, the solution is far from the optimal solution. Inverse con- 
dition number of the eigenvector matrix is 1.26 x 10 -5 indicating the presence of degenerate 
modes. 

This degeneracy is confirmed when the optimization was run using the robust form for 
the gradient computation. The optimized cost function now has a value of 2.38x 10 -12 which 
is much smaller than that achieved under the diagonalization algorithm. In addition, inverse 
condition number of the eigenvector matrix is now 9 x 10 -6 indicating that degeneracies 
occur through the design optimization. 


3.11 Conclusions 

Numerical algorithms for computing matrix exponentials and integrals of matrix exponen- 
tials have been developed to handle cases where the system matrix has defective eigenval- 
ues. Special formulation of these algorithms enables reliable and efficient computation of X 
and M. These algorithms have been incorporated into the computer-aided design package 
SANDY for synthesizing optimal output -feedback controllers. Numerical optimization com- 
bined with the given algorithms in the evaluation of the cost function and its gradients with 
respect to the controller design parameters is shown to have well-behaved convergence even 
when the closed-loop system becomes degenerate. Rehability of the algorithm is further 
demonstrated using typical design problems encountered in control of flexible structures. 
Clearly, this algorithm, when combined with a previous one based on diagonalization, would 
enhance significantly the overall reliability of optimal control synthesis using parameter op- 
timization, thereby providing an effective automated design environment for multivariable 
control synthesis. 

To further improve the efficiency of the robust algorithm, a hybrid form combining the 
benefits of both the diagonalizat ion-based form and the robust form has been developed. 
There are two major components to this algorithm. The first one is a means of reliably 
decomposing the system matrix into two parts: a diagonal block containing nondefective 
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modes and diagonal subblocks containing the defective portions. Details are presented 
in Appendix A. The actual calculation of X and M are performed in the second step, 
consisting of elements in the diagonalization-based form, the robust form, and new cross- 
term calculations. Details are discussed in Appendix B. 
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Chapter 4 

Control Design for a Large Space 
Structure 


4.1 Introduction 

In this chapter, we present a control design problem that poses significant difficulties to 
the early design algorithm based on diagonalization. The system model contains numerous 
densely packed flexible modes primarily due to the symmetry in the configuration of the 
structural design. There is a high likelihood of numerical degeneracies due to round-off in 
a finite-precision machine. The associated eigenvector matrix tends to be ill-conditionned 
and this problem occurs frequently during the design optimization. 


4.2 JPL Large Space Structure 


Defective degeneracy in low-order systems can be easily identified and the problem possibly 
remedied through simple changes in the model parameters. For example, transverse Dryden 
turbulence filters contain a pair of degenerate real roots; these filter poles can be perturbed 
slightly to produce a non-defective system with virtually no loss of accuracy in the calcu- 
lation of power spectral densities. Most linear time-invariant controller state matrices can 
be represented in a canonical form with 2x2 blocks along its diagonal of the form 


0 1 
— cj 2 — 


(4.1) 


Clearly a choice of £ = 1 would lead to a degenerate subblock. 

In control problems for large flexible mechanical systems such as space structures, causes 
of eigenvalue degeneracies are usually more subtle in nature than the simple case presented 
in Section 3.9 for a two-mass-spring system. The JPL large space structure has been 
carefully designed to simulate a lightweight, non-rigid and lightly damped structure in a 
weightless environment [28]. The structure itself resembles a large antenna with a central 
boom-dish apparatus and an extended dish consisted of hoop wires and 12 ribs. There are 
two torque actuators (labelled HA1 and HAIQ) on the boom and dish structure to control 
the two angular degrees of freedom in pointing maneuvers, and force actuators at four rib 
root locations (labelled HAl, RA4, RA7 and RAIO) for vibration control. From the point of 
view of control design, it is a challenging problem, since the plant has many closely spaced 
modes and is of reasonably high order. There are a total of 30 modes in the basic structural 
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Figure 4.1: Antenna Structure 
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model. The flexible modes are lightly damped with damping ratios ranging from 0.007 to 
0.01. The two rigid-body modes have a damping ratio of 0.12. Our design concept is to use 
two available angular displacement sensors HS 1 and HS 10 of the boom-dish apparatus and 
the two torquers HAl and 77.410 collocated with these sensors for control synthesis. With 


Table 4.1: Open Loop Modes of the Antenna Structure 


Eigenvalue 

Damping 

Freq (Hz) 

-0.09500 

± 

0.7860* 

0.120 

0.12600 

-0.08575 

± 

0.7093* 

0.120 

0.13704 

-0.02802 

± 

4.0024* 

0.007 

0.63701 

-0.02929 

± 

4.1844* 

0.007 

0.66598 

-0.07405 

± 

10.583* 

0.007 

1.68434 

-0.07405 

± 

10.583* 

0.007 

1.68434 

-0.11310 

± 

10.616* 

0.007 

2.57123 

-0.11785 

± 

16.384* 

0.007 

2.67929 

-0.21365 

± 

30.520* 

0.007 

4.85749 

-0.21365 

± 

30.520* 

0.007 

4.85749 


this selection, 20 of the flexible modes associated primarily with the rib motion become 
uncontrollable and unobservable. These modes are removed by modal truncation from our 
plant synthesis model. Eigenvalues of the remaining 10 modes are shown in Table 4.1. 


4.3 Controller Design 


An optimal low-order controller is designed to dampen vibration of the antenna in response 
to external excitations. To evaluate the effectiveness of the control system, we perform 
the following test. The entire structure is agitated using the two boom-dish actuators for 
the first 6.4 seconds with an applied torque in the form of a square wave of 0.8 second in 
width and with an amplitude of 1 N-m. The control system is then activated right after 
the excitation has been removed, and responses of the excited structure at the sensors are 
examined. The design objective is to damp out the induced vibration as fast as possible 
without excessive use of controls. Note that the natural responses of the structure will take 
about a few minutes to decay to zero (Figures 4.2 and 4.3). 

For practical implementation, the controller design is choosen to be of 6 th order. The 
controller structure consists of a pair of tightly coupled second order systems each with 


dynamics of the form 


0 1 
— u ) 2 


a pair of actuator lag states, 
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The two first-order lag states in the controller model serve as roll-off filters, limiting the 
control bandwidth to less than 50 rad/sec. In the design optimization, we have a total of 
28 design variables: 16 in the controller A matrix and 12 in the B matrix. The objective 
function for design optimization consists of a sum of weighted // 2 -norms of physical response 
variables observed at different locations of the structure. It is of the form 

if 12 2 1 

J(t f ) = Lira - 2 VEQ iEa [yKt S )\ + £ RjE a [«?(*/)] j (4.3) 

Note that the expectation operator E a [*] is for a system destabilized by a factor a. Table 4.2 
lists the design variables j/i and their corresponding penalty weightings Q t . Also given in 


T able 4.2: Design Variables: JPL Antenna Structu re 


Variable 

Qi 

Description 

RS 1 

4100 

Rib #1 root velocity 

i R'S4 

3950 

Rib #4 root velocity 

RS 7 

3975 

Rib #7 root velocity 

RS 10 

4050 

Rib #10 root velocity 

HSl 

16500 

Hub angular velocity 

HS 10 

15600 

Hub angular velocity 

RSI 

1100 

Rib #1 root displacement 

RS4 

1050 

Rib #4 root displacement 

RS7 

1150 

Rib #7 root displacement 

RS10 

1025 

Rib #10 root displacement 

HSl 

3900 

Hub angular displacement 

HS10 

4100 

Hub angular displacement 

Variable 

R, 

Description 

HAl 

41 

Hub torque actuator 

HA10 

40 

Hub torque actuator 


the table are the control design weightings Rj for the actuators HAl and HA10. Responses 
in the above objective function are evaluated to random disturbances of unit white-noise 
spectra applied simultaneously at all the hub and rib actuators. One may notice that the 
values of related weighting terms are perturbed about a nominal. With this nominal value 
as the weighting factor, the perturbation was for avoiding degenerate modes in the optimum 
design. 

The design optimization begins with the following arbitrary initial guess on the controller 
matrices A and B, 


' -50 

0 

1 

0 

0 

0 ' 


‘ 0.1 

0 ' 

0 

-50 

0 

0 

1 
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0 

0.1 

0 

0 

0 
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; B = 

0 

0 

0 

0 

-2 

-1 

0 

0 

0 

1 

0 

0 

0 

0 

0 
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0 

0 

0 

0 

0 

0 

-4 
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1 

0 . 


A destabilization factor a of 0.071 was used to ensure that all the closed-loop eigenvalues 
have a real part less than -0.071. The optimization fails to converge when a destabilization 
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factor of greater than 0.075 was selected. This difficulty seems to be in moving the modes 
at 1.68Hz under this controller configuration, implying that additional degrees of freedom 
must be added to the controller structure given in equation (4.2), in order to improve the 
performance further. 

While the optimization convergence itself took 13.5 hours on a VAX/VMS Workstation 
DEC-3500, the proposed algorithm for the calculation of the objective function and its 
gradients with respect to the design parameters is robust and leads to well-behaved design 
convergence. The final optimal values of the A and B matrices are 


50 

0 

2.874 

-2.270 

1.7322 x 10“ 3 - 

2.2131 x 10" 4 

0 

-50 

1.225 

0.7825 

6.551 

-1.037 

0 

0 

0 

1 

0 

0 

0 

0 - 

15.73 

-0.8799 

0 

0 

0 

0 

1.560 

0.2256 

0 

1 

0 

0 

2.400 

-1.269 

-13.62 

-0.9810 



r 

5.343 

-1.2310 x 10" 4 ' 




6.2118 x 10" 4 

4.783 



B = 


2.701 

-8.1595 x 10~ 4 




2.221 

9.3152 x 10“ 4 





-0.5147 

5.379 





1.614 

-1.252 



Closed-loop responses of the sensor and control variables corresponding to this design are 


Table 4.3: Boom-Dish-Controller Closed Loop Modes 


Eigenvalue 

Damping 

Freq (Hz) 

-0.086899 

± 

0.6588i 

0.1308 

0.1058 

-0.089071 

± 

0.7410i 

0.1193 

0.1188 

-0.3165 

± 

3.624i 

0.0870 

0.5790 

-0.2528 

± 

3.790i 

0.0666 

0.6045 

-0.2162 

± 

4. 1 12i 

0.0525 

0.6553 

-0.2056 

± 

4.1851 

0.0491 

0.6669 

-0.074193 

± 

10.58i 

0.0070 

1.684 

-0.074589 

± 

10.581 

0.0070 

1.684 

-0.1168 

± 

16.1 5i 

0.0072 

2.570 

-0.1253 

± 

16.83i 

0.0074 

2.678 

-0.2142 

± 

30.52i 

0.0070 

4.857 

-0.2143 

± 

30 . 52i 

0.0070 

4.857 

-49.99 



1.000 

7.956 

-49.99 



1.000 

7.956 


shown in Figures 4.2 and 4.3. The controlled responses decay to zero in about 20 sec after 
the excitation has been removed. Notice that the control torques are within the desired 
limits of 1 N-m; the results are obtained through the control design weights Rj in Table 4.2. 
This design example demonstrates the usefulness of a design algorithm for robust low-order 
controllers using parameter optimization, and the accompanying improvement of solution 
reliability using the algorithms described in Sections 3.6 and 3.7 for degenerate systems. 
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Figure 4.2: Hub Responses to Structure Excitation 
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Figure 4.3: Hub Responses to Structure Excitation (cont.) 
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Chapter 5 


Closed-Loop Transfer Recovery 


5.1 Introduction 

Multivariable control has evolved over the last decade. Methods for synthesizing control- 
laws using optimal control from linear quadratic regulator to linear quadratic gaussian de- 
signs have been studied extensively. Robustness issues of LQG have led to the development 
of the concept of loop transfer recovery [1]. Recent work by Saberi et. al has fundamen- 
tally examined the loop transfer recovery based on a system theoretic approach [11]. In 
this chapter, we review briefly the concept of closed-loop transfer recovery which is an ex- 
tension of the loop transfer recovery method for the recovery of closed-loop input/output 
behaviour. Observer-based controller designs for CLTR will be presented to illustrate the 
design concept. A new alternate approach based on numerical optimization is developed 
to solve CLTR designs for low-order compensators. This approach is then applied to the 
synthesis of a high-bandwidth control system for a rotorcraft in Chapter 7. 


5.2 Analysis of Closed-Loop Transfer Recovery 

For performance and robustness, multivariable control synthesis usually begins with the 
design of a state-feedback controller. A state-feedback design based on LQR techniques 
involves only the solution of Riccati equations and provides a convenient and efficient way to 
examine the control problems in terms of achievable performance and robustness. However, 
state-feedback design is not practical due to its requirement of noise-free measurement of 
all the states, but such a design can provide a target response for output-feedback designs 
using loop transfer recovery (LTR). Recovery of closed-loop transfer functions achieved 
under state feedback is classified under the category of closed-loop transfer recovery. This 
concept has been developed in [11]-[12] using full and reduced-order Luenberger observers. 
Designs can be solved using Riccati equations. This approach is later extended to arbitrary 
low-order compensators using numerical optimization. 

The premise of CLTR is simple. Consider a linear time-invariant plant 


x(t) 

= Fx{t ) 

+ 

r™(t) 

+ 

Gu (t) 

(Dynamics) 

z(t) 

= H c x{t) 

+ 

Ffytjj VJ ( t ) 

+ 

DcuU(t) 

(Criteria) 

y(t ) 

- H s x(t ) 

+ 

DswWit) 

+ 

DguuiV) 

(Sensors) 


where w(t) represents the command/disturbance input vector and u(t) is the control input 
vector. We assume that a state-feedback design u = Kx has been obtained that achieves 
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the desired closed-loop properties in the transfer function T zw (s) between the input w(t) 
and the output z(t). Many different design methods can be used, e.g. H 2 - optimal control 
[15], //“-optimization [4], and others [33]. A state-feedback design is totally characterized 
by the gain matrix K and delineates the achievable performance. The actual performance is 
often limited due to a restricted set of measurements and the associated sensor noises. The 
problem reduces to finding output- feedback controllers to recover as much as possible the 
performance achieved under state feedback. The solution can be found using the concept 
of closed-loop transfer recovery with the recovery error defined by 

E(s) = T° w {s) - T zw (s ) = T zu (s)M(s ) (5.1) 

where T zw (s ) and T° w (s) axe the closed-loop transfer functions between the disturbance 
input w(t ) and the criterion output z(t) of a state-feedback and observer-based output- 
feedback designs respectively. The transfer function T zu (s) is simply the closed-loop transfer 
function between the control input u(t) and the criterion output z(t ) under state feedback 
(independent of the output-feedback design). The transfer function M(s) takes on different 
forms depending on whether a full or a reduced-order Luenberger observer has been used 
[11]. The structure of the auxiliary transfer function M(s) is defined explicitly in Sections 5.5 
and 5.6. The CLTR procedure for observer-based controllers addresses the minimization 
of the norm of ||M(s)|| whose solutions can be found from Riccati-based methods in either 
H 2 - or //“ based optimization. However, minimization of the norm of the actual error 
function ||/?(s)|| can only be done using numerical optimization schemes to be discussed in 
Section 6.2. 


5.3 The Special Coordinate Basis 

Analysis of closed-loop transfer recovery can be effectively performed when the system model 
is transformed into a Special Coordinate Basis (SCB) [10]. Under the SCB transformation, 
many system properties can be identified: right and left invertibility, system invariant zeros, 
infinite zero order structure and the related geometric subspaces. 

Appendix C presents an alternative approach underlying the conceptual development of 
the SCB transformation following similar notation as given in [10]. Other useful applications 
of the SCB properties are found in H 2 and H°° optimization [33]. 

Let’s consider for simplicity a strictly proper system. The basic notion of the S.C.B. is 
to divide up the system states into internal states and output states. The output states are 
either themselves outputs or are inputs to a cascade of integrators leading to the outputs. 
The state at the top of each cascade will either have one control term entering into its 
state equation, or it will not have any control input at all. Thus each cascade of output 
states will belong to one of two classes. The internal states are the remaining system states 
that do not fall into this category. The internal states can be collectively thought of as 
providing feedback to the system from the outputs (Figure 5.1). The class of internal states 
is subdivided into those states that are controllable with those control inputs not present 
in the dynamics of the output states, and those states which are not controllable with these 
inputs. This partition is slightly analogous to that of the output states. Note that if an 
“internal” state has an input term that directly affects an output state, then this input can 
simply be replaced by the corresponding output state and hence leaving this internal state 
with no direct input term. In summary, 
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Internal States 

Output States 

Xa • 

No direct input 

x c : 

Direct input v 
{u and v disjoint) 

VS ■ 

Direct input u 

Vb ■■ 

No direct input 


Invariant zeros are the eigenvalues of the partition of the SCB transformed system matrix 
corresponding to the states x a . Left and right invertibility are easily visualized in the SCB 
transformed space and are not to be confused with controllability and observability. Left 
invertibility can be defined as follows: for a given set of outputs, there exists a unique 
set of inputs that generates these outputs. A system is left invertible if the states x c are 
non-existent. Similarly, it is right invertible if the states y b are non-existent. 

The condition of left invertibility is not the same as that of observability. All the output 
states y/ and y b are observable, so a lack of observability may only exist in the internal 
states x a and x c . By the same token, some of the output states y b and/or some of the 
internal states x a may be uncontrollable. 

The SCB transformation allows designers to examine the general system input-output 
properties for a given set of inputs and outputs. It should be noted that the inputs can 
be either controls and/or disturbance inputs, and the outputs can be either sensor and/or 
criterion outputs. 

In the context of closed-loop transfer recovery, conditions for either exact or asymptotic 
recovery depend on the system input-output properties between the disturbances w(t) and 
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the sensor outputs y s (t) for the state-feedback stabilized system 

v' . f x(t) = (F — GK)x{t) + Gu{t) + Tw(t) , . 

V ‘ W ' \ Vs{t) = H s x(t ) + D m u(t) + D sw w(t ) ^ ' 

This analysis applies only to recovery designs with a full or reduced-order observer-based 
controller. Issues related to controllability are implicitly resolved with the existence of 
a suitable state-feedback law. In the next section, we recall some basic conditions for 
recoverability obtained in [11] for observer-based controllers. 

5.4 Conditions for Exact and Asymptotic Closed- Loop Trans- 
fer Recovery 

As previously mentioned, we need to examine the SCB properties of the system E VsW be- 
tween the disturbance input w(t) and the sensor output y s (t). It has been shown [11] that 
recovery is possible if this system T, VaW is left invertible and has stable invariant zeros. For 
the case of full-order observer-based controllers, exact recovery is possible if furthermore 
the system E yaUI under the SCB transformation has no yj states; otherwise, one can only 
achieve asymptotic recovery. In another word, full-state feedback design is not recoverable if 
the system E ysU) has x c states (i.e., not left- invertible) or nonminimum phase (i.e., invariant 
zeros on the right-half plane). 

In the case of reduced-order observer-based controllers (i.e., Luenberger observer), we 
have the same recovery conditions as in the case of full-order observer-based controllers. 
However, exact recovery is still possible even when there are output states y/. The exact 
recovery condition requires simply that each output state in yj is separated from the inputs 
u{t) by at most one integrator. If more than one integrator exists between them, then the 
recovery can only be achieved asymptotically. 

A more complete discussion on recoverability of observer-based controllers can be found 
in [11], For non-recoverable systems, the nonzero values of the auxiliary transfer function 
output M{s) are transformed to the actual error 

E(s)=T zu (s)M(s) 

which could be quite large for a nonzero M(s). For the non-recoverable case, the CLTR 
design problem is more appropriately addressed through minimization of the actual recovery 
error E(s), which makes the use of numerical optimization necessary to determine the design 
solution. This design approach is presented in Chapter 6. 


5.5 CLTR Design with a Full-Order Observer-Based Con- 
troller 

Consider the recovery error achieved under a full-order observer-based controller [11], 

E f (s) = T<J>(s) - T zw (s) = T zu (s)Mf(s) 

with T zu (s) is the closed-loop transfer function between the control u(s) and the criterion 
output z(s) under the state-feedback design. Similarly, T zw (s) is the closed-loop transfer 
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between the disturbances w(s) and the criterion output z(s). The transfer function T<f > (s) 
corresponds to the closed-loop system with a full-order observer-based controller. 

Consider the following plant and full-order observer-based controller, 


Plant 

x(t) = Fx(t) +Tw(t) + Gu(t) 

z(t ) = H c x{t) + D^w^) + £>cuu(4) 

y s (t) — HsXjt) + DawW(t) + Dstjujf) 

Full- Order Observer-Based Controller 
&{t) = (F - K ot>3 H s )x{t) + Gu{l) + K^yjt) 

u(t ) = —Kx(t) 


where K is the state-feedback gain matrix and K 063 is the observer gain matrix to be 
determined under the CLTR procedure. The auxiliary error Mj(s) can be described directly 
in terms of the following state model 

M f (s ) = K(sl -F + K^HsY'iT-K^Dw) 

Mj(s) = (r T - Dj w K obsT )(sI -F T + HfK obsT )- 1 K T 

The system transfer function Mj(s) can be described by the following auxiliary system 

j i{t) = F T x(t) + Hju(i) + K T w(t) 

\ z(t) = r T x(t) + Dj w u(t) 


k T T 

under the state-feedback control u(t) = —K°° s x(t). Synthesis of the observer gain K°° s 
can be performed through the minimization of the H 2 - norm of the criterion output z(t) as 
in an LQR state-feedback problem. Namely, 


mm 

K obs 


roo roo j- 

J z T (t)z(t) dt = min J |x T (f) u T {t ) J 


rr 1 


TDl 


d r T n d t 

l ^SW 1 Ly SW J - y $‘U) 


x(t) 

u(t) 


dt 


This problem is a singular // 2 -optirnal control due to the fact that the weighting matrix 
D sw Djw is usually singular. One can solve this problem by introducing a small perturbation 
into the singular portion of this matrix. To do this, we first perform a singular value 
decomposition on D sw as follows, 

D sw = U d ZdVE =» = UdY.dY.IuI 

0 0 ••• ] 

a\ 0 ••• 


0 

Singularity of a matrix is usually defined relative to its maximum singular value as in the 
definition of the condition number. Let’s choose a small perturbation equal to with e in 


wneie 




Y d Y t d = 


0 

0 
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the vicinity of the machine precision. If this number is greater than < 7 ^ for some m, then 
one would perturb the matrix D sw Dj w by 


R — [D.sw D SVJ }perturbcd — Dsw E sw ” 1 " €(T\Uq 


[0]( m -l)x(m-l) 0 

0 [I 


ul 


This problem can be solved as a regular LQR problem and solution obtained from an 
algebraic Riccati equation 

0 = 5 [f t - HjR'- l D sw R T ) + (f~ TDj w R'- l H a ) S + T (/ - D SW R'-' D sw T t ) F t 

—shJr'~ 1 h s s 

with K obs = R'~ l ( H s S 4- D sw T T y Let’s define E(s,e) = T£/ > (s,e) - T zw (s ) associated 

with the family of closed-loop transfer functions T</ > (s,e). One can examine the asymp- 
totic properties as e — > 0. If the conditions for exact recovery are satisfied, then for e — ► 0, 
the gain matrix K 063 will converge to a finite gain matrix as E(s, e) — > 0. Otherwise, some 
of the gains in K ohs will approach infinity. Closed-loop responses of the state-feedback 
design associated with infinite zeros of the system E yilW will not be recovered asymptoti- 
cally while those corresponding to the finite stable invariant zeros will be exactly recovered. 
In the case of asymptotic recovery, although E(s, e) — > 0 as e — » 0, one must pick an e 
representing a desired level of recovery along with a reasonable finite gain for K^ 3 . 


5.6 CLTR with a Luenberger Observer 

In the Luenberger observer design, some of the measurement outputs are assumed to be 
noise free. For a given state-feedback matrix K, a controller based on the Luenberger 
observer is of the form, 

v(t) = Lv(t) + Miy s (t ) + M 2 u(t) 
x(t) = Pv(t) + Jy s (t ) 
u(t) = —Kx(t) 


with 

QF -LQ = M\H S , M 2 = QG, 

and 

JH S + PQ = I n 

Note that n is the order of the plant model. 

Without loss of generality, we can always transform the plant model such that the 
noise-free outputs of y s {t ) become a subset of the system states (see [11], Section 4.2). 
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Clearly, an observer is needed to estimate the states X 2 H). Defining an observer gain matrix 
K° bs = [K$ s , K£j s J whose partitions correspond to the measurement outputs yo(t) and 
y\ (f) respectively, we obtain a controller with a reduced-order Luenberger observer of the 
form [12], 


m 

x(t) 

u(t) 


Aorv(t) + (g 2 - K$ S G ,) u(t) + K$ s y 0 (t) 
+ (f 2 1 - K°\ S F\\ + AorKfr) yi(t) 


0 

In — p+THo 

—Kx(t) 


v(t) + 


Ip-mo 

J£ObS 


y i(0 


where 

Aar = F22 - K% S H 3 ,02 ~ K$ S F 12 

We also partition the state feedback gain matrix K = [K\, K 2 \ in accordance with the state 
partitions X\ (t) and x 2 (t). The recovery error can then be written as in [11], 


E r (s) = T<J>(s) - T zw (s) = T zu (s)M r (s) 


where T^ > (s) is a closed-loop transfer function incorporating a reduced-order observer- 
based controller, and 

M r {s) = K 2 (si - A r + K^ay 1 (B r - K? s D r ) 

or 

Mj ( s) = {bJ - Dj K? s y {si -Aj + cj K? s y ~ l Kj 

where 

At — F 22 , B t — T2, 

Bswft 

r 1 

As in the case of full-order observer-based controller, an auxiliary system can be constructed 
such that its closed-loop transfer function under a “state-feedback” law is equal to Mj(s). 
It is given by 

x r (t) — Ajx r (t) + Cj u(t) 4- Kjw(t ) 
z(t) = Bj'xr(t) + Dju{t) 

and u(f) = —K^Xrit). 

The reduced-order observer design gain matrix Kf* 3 can be obtained from the mini- 
mization of the // 2 -nomn of the criterion output z(t). If the associated LQR problem is 
singular, one can perturb the weighting matrix D r Dj in a design approach similar to the 
that presented in Section 5.5. As in the full-order observer case, closed-loop responses 
corresponding to finite stable invariant zeros will be exactly recovered. Furthermore, re- 
sponses associated with first-order infinite zeros will also be exactly recovered. Responses 
corresponding to infinite zeros of second and higher order will be asymptotically recovered. 
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5.7 Conclusions 


In this chapter we introduced the concept of closed-loop transfer recovery. Analysis of 
recoverability using observer-based controllers can be done using system properties derived 
from the SCB transformation. A procedure based on Riccati solutions is given for the design 
of full and reduced-order observer gain matrices. 
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Chapter 6 


CLTR Using Numerical 
Optimization 


6.1 Introduction 

Although numerical optimization provides a powerful design tool for the synthesis of multi- 
variable control systems, its usage is still limited. Difficulties often encountered in optimal 
control synthesis are related to how one formulates the design problem to account for all 
the basic requirements in stability, performance and robustness. Even in // 2 ///°°-optimal 
control, extensive design trade-offs are often needed before a satisfactory optimum design 
can be found. Some of the trade-offs reside in the selection of a quadratic performance in- 
dex. Here, one needs to identify, through successive trials, the appropriate design criterion 
variables, the control variables, and their relative effects on the overall design. The process 
tends to be computationally intensive and time consuming if one considers direct synthesis 
of low-order controllers using numerical optimization. However, if the control problem can 
first be solved under the setting of state feedback whose solutions can easily be obtained 
from the algebraic Riccati equation, one can then proceed to the design of output-feedback 
controllers using the concept of closed-loop transfer recovery. The merits of this approach 
are that few design iterations are needed in obtaining an output-feedback controller, espe- 
cially when a low-order controller is derived from the order reduction of an observer-based 
controller that achieves exact recovery. 

In this chapter, we address the problem of closed-loop transfer recovery using an ar- 
bitrary controller structure. As indicated in the previous chapter, analysis and design of 
closed-loop transfer recovery are well developed in the setting of observer-based controllers 
[11]-[12]. These results form the baseline designs to CLTR using numerical optimization. 
Performance and robustness achieved under observer-based controllers serve as guidelines 
in defining a satisfactory output-feedback design of lower order. A tentative low-order con- 
troller design can be obtained from order reduction of an observer-based controller using a 
balanced truncation method. Hankel singular values of the observer-based controller model 
can be used to determine the possible order of the dynamic compensator in numerical op- 
timization. The advantage of numerical optimization is that it allows designers to seek 
controllers of a specific (practical) order and structure in the recovery of the closed-loop 
characteristics achieved under state feedback. The numerical scheme can further be used to 
re-optimize the recovery error over a given control/command bandwidth using frequency- 
shaped filters. 
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6.2 Problem Formulation 


In this section, we define the problem associated with closed-loop recovery using an arbitrary 
dynamic compensator. Under this setting, there is no longer the simple relation for the 
recovery error as defined in equation 5.1. 

For a dynamic compensator of the form 


i . f x c (t) — A c x c {t) + B c ys{t) 
' c ' ( u(t) = C c x c (t ) + D c y s (t) 


analytical results for closed-loop transfer recovery analysis and design are not possible 
in general. Closed-loop transfer recovery problems presented in [11] and [12] apply only 
to observer-based controllers. Given a recoverable system for an observer-based design, 
a controller with the general structure shown above is also likely to recover the desired 
state-feedback properties. 

With H 2 optimal control, H°° optimization, pole placement, or other design techniques, 
one can determine a set of state-feedback gains K that meet desired closed-loop character- 
istics through the control u(t) — —Kx(t). The resulting closed-loop system 


f x f (t) = (F -GK)x f (t)+Tw{t) 

[ z(t) = (H c - DcuKjxftt) + Dcwwit) 


defines the ideal transfer function T zw (s) to be recovered with an output-feedback design. 

In the design approach based on direct optimization for CLTR, the quantity being 
minimized is the actual error 


E(s) = T zw {s) - T° zw (s) 


rather than the auxiliary transfer function Mj(s ) (equation 5.1). Typically the disturbance 
input w(t) is modelled as white noise. Frequency shaping can be introduced into w(t) using 
additional filter dynamics 


f x w (t) = F w x w (t) + T w r](t) 

( w(t) = H w x w (t) + D w rj(t) 


where rj(t) are white noises. The overall synthesis model, including the noise shaping 
filters, the state-feedback target system dynamics, the open-loop plant, and the controller 
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dynamics is given by 
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where e(t) is the recovery error, and y e {t) is the measurement output vector augmented 
with the states of a dynamic output-feedback compensator of order r. Feedback control is 
conveniently implemented in a static output feedback form 


u(t) 

U c (t) 


— Co2/e(f) — 


D c C c 

B c Ac 


Ve{t). 


The CLTR problem for a dynamic output-feedback compensator reduces to finding a gain 
matrix C 0 that minimizes a certain norm of the recovery error, or 


min JcUr(C 0 , t f ). 

ts o 

For the H 2 norm of the recovery error, the disturbance inputs are treated as white noises 
with E[rj(t)r](T) T ] = W a S(t — r). The corresponding H 2 cost is given by 

Jdtr = hm^E [ e T (tf)e(tfj\ . 

By Parseval’s theorem, it has the equivalent form 






2 


du 


in the frequency domain. It can be extended to include a frequency weighting W(ju) as 
follows 

1 A°° 

Jdtr = -J Q \\[T zw (ju) - T% w (juj)} W(ju) II 2 du>. 

Due to the difficulty in finding an analytic solution for a general controller structure, a 
minimum J c Ut is often obtained by nonlinear optimization techniques. 
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Difficulties in numerical optimization for the CLTE approach are associated with the 
following two factors. The first is that poles of the closed-loop system in the output- 
feedback design are optimized toward those of the “reference” state-feedback system. While 
the eigenvalues are coming closer together, the modes themselves remain distinct because the 
dynamics of the state-feedback and the output-feedback systems contained in the synthesis 
model remain distinct. The typically large size of the overall system matrix brings in the 
second factor. Roundoff in computations involving large system matrices could cause non- 
degenerate modes with the same eigenvalue to appear degenerate. With degenerate modes, 
or correspondingly, defective eigenvalues in the system matrix, algorithms for computing 
the cost function and its gradient based on an eigenvalue-eigenvector decomposition will 
break down. Since a degenerate system matrix is likely in the CLTR procedure for most 
practical problems, one must employ techniques for computing the cost and its gradient 
that are immune to this condition. 


6.3 Feedforward Controller in the CLTR Design Procedure 

We present a design approach for feedforward gains in numerical CLTR. Determination of 
the feedforward gains in the optimization process would certainly provide the best responses. 
However, because the feedforward gains do not affect stability, one need not include them in 
the controller optimization process, thus reducing the number of parameters involved and 
saving computation time. 

A particular feedforward structure is chosen here for convenience. Given a feedback 
controller design, the basic idea is to construct a feedforward gain to achieve the desired 
steady-state relations between command inputs and the commanded outputs. Although 
we limit the scope to constant gain matrices, the resulting solutions are often non-unique. 
First and foremost, we review the computation of feedforward gains for the state feedback 
case. In the subsequent sections, we consider observer-based controllers followed by the 
case for a general controller structure. 


6.3.1 A Feedforward Controller under State— Feedback 


This is the simplest case of feedforward controller design. We write the system dynamics 
and commanded output in a single set of equations 


x(t) 


ycmd (0 



(F - GK) 
{Cand D cmc iK') 
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Dcmd 


x{t) 

^•cmdify 


( 6 . 1 ) 


Here «cmd(i) is the direct control input derived only from the command feedforward gains 
acting upon y re f(t)- A feedforward controller is then determined such that the system 
outputs ycmd{t ) match the command reference inputs y re f(t) in steady-state, i.e., when 
x(t) = 0. Solving equation (6.1) for steady-state values of the states and the control inputs, 
we have 
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Figure 6.1: State Feedback System with a Feedforward Gain 


where x ss is the steady-state value of x , and 

L 1 _ I" ( F-GK ) G 1 _1 [ 0 ' 

Af {Gcmd G cm dK') Ocmd I 

From this, we can deduce that u cm d(t ) = Ny re f(t ) and therefore Kjf = N. Very often the 
matrix 

(F - GK) G 
{Gcmd D cmc [ 

is not full rank, therefore we must take the pseudoinverse rather than the actual inverse. 
This signals the existence of more than one solution to the feedforward gains (we have more 
degrees of freedom than necessary). This feature comes in use when we pick convenient 
solutions to some of the cases that follow. 

Considerable simplification can happen when one considers command inputs to integral 
states. There are two similar, but distinct cases involving integral states recognized here. 
In the first case, the integral states are part of the plant model 


' *(t) 1 [ F 0 

x T {t) j F 1 0 




Assume there is no direct control input into either the integral states or the commanded 
output ycmd (G r = 0 and Dcmd = 0). Partitioning the state feedback gain matrix into 
\K K r ], the overall system would look like 

x(t) ] F-GK -GK 1 G 1 [ x(t) 

x\t) = F 1 0 0 x ! {t) . (6.3) 

Vcmdii) 0 Gcmd 0 ^cmdif) 

Again, the steady-state solution would correspond to x = 0, x 1 = 0, and y cm d = Vref ■ Let’s 
consider an alternative solution u CTn d = K 1 C^ d y Te j or Kjf = K T C^ d . (It is assumed, 
without loss of generality, that Ccmd is full rank). With this feedforward gain, we have 
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x ss = 0, and x[ s = C^yre/. Other nontrivial solutions (nonzero x 3S ) of equation (6.3) are 
possible. 

The second case involves plant dynamics that are augmented with the integral states 
i / (t) = y CTnd (f) —y re f(t). The combined equations for the augmented system dynamics and 
the commanded output are given by 
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In steady-state, we have 
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While this can be solved for Kff as before for Ucmd = Kffy re f } it must be noted that the 
command path through the integral states will govern the long-term command response. 
The structure of the above equations guarantees the existance of more than one nontrivial 
solution (x 3S ^ 0) for Kff when such a solution exists, allowing a choice of transient 
responses to command inputs. Often Kff is set to zero. 

6.3.2 A Feedforward Controller with Observer— Based Controllers 

Design of the feedforward gain matrices Kff and G c (as shown in Figure 6.2) will differ 
slightly between a feedback controller design based on an observer-based controller structure 
and one with an arbitrary dynamic compensator structure. 

A direct feedforward input into the controller dynamics provides an “anticipation” factor 
to the command inputs. An observer-based controller using the CLTR procedure possesses 
an inherent feedforward controller structure. This is due primarily to the fact that the 
observer has a direct input term from the control u(t). Let’s examine in more detail the 
influence of a command input in an observer-based controller design. The observer-based 
controller is given by 


x c (t) = A c x c (t) + B c y s (t) + G 2 u(t) 
x(t) = Px c (t) + Jy 3 (t) 

where the control output is 

u(t) = Kx(t) + 

'LL cm d{t) 

— KPx c (t) + KJy s (t) + Uamdit) 

— C c X c {t) “1" G c y s (t) + licmd(^) 

— C c X c {t) + D c [H s x(t) + DguU^t)] + Ucmd{t) 
= C c Xc(t) + D c H s x(t ) + (I + D c D su )u C md{t). 
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Figure 6.2: Closed-Loop System with a Feedback/Feedforward Controller 


where we deine C c = KP and D c — KJ. For smplicity, we assume D C D SU = 0. The 
combined closed-loop dynamics and command ouput equations are given by 
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In steady state, we have 
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We have Ucmd{t) = Ny re f{t ) or K/j = N. Within the context of Figure 6.2, it is straight- 
forward to see that G c — G 2 Kff. 

Like the state-feedback case, some simplification is possible if one considers problems 
with integral control 
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Note that the control is given by u(t) = C c x c (t) + D c y s (t ) + Df.x r (t). This would result in 
a closed-loop system of the form 
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As in equation (6.4) for the state feedback case, the above form allows more than one 
solution to the command feedforward u cm d = Kjfy re /, with K/f = 0 often being used. As 
before, only the short-term response is governed by this feedforward. 


6.3.3 A Feedforward Controller with Non-Observer— Based Controllers 

We can now proceed to the design of the feedforward controller for an arbitrary feedback 
compensator of the form, 

x c (t) = AcX c (t ) + B c y„(t) 
u(t ) = C c x c (t ) -I- D c y s {t ) 

The overall closed-loop structure is exactly as shown in Figure 6.2. Because there is no pre- 
set distribution from the control Ucmd to the controller, determination of the feedforward 
gains K/f and G c will have more degrees of freedom due to the nominal independence of 
the two matrices. 

We will cover three common approaches for designing the feedforward gain matrices K // 
and G c ■ In the most general case, one simultaneously creates both matrices. We first write 
state and output equations for the overall closed-loop system 
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where Ucmd(£) = G c y re j(t ) and we again assume D C D SU = 0. The above system is not 
square, with more degrees of freedom than equations. Thus, one can solve the above using 
a singular value decomposition on the system matrix 


U f Y, f V$ = 
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by constructing a pseudoinverse. The form of the pseudoinverse is V^Ept/J, where E^ is 
generated by transposing the singular value matrix Ep and replacing each nonzero singular 
value with its reciprocal. Applying this to the steady-state response, 
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The other two methods depend on pre-selecting either of the matrices K/f or G c . In 
direct optimization of a controller, one often has the commanded outputs coinciding with 
several of the sensed outputs, in which case it would be reasonable to make —G c equal to 
selected columns from B c . The overall system becomes 
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The pre-selection of K jj happens in the context of CLTR. Since one is trying to recover 
state-feedback characteristics, the steady-state controller response should be like that from 
the state-feedback gains. Therefore it is appropriate to use the K ff from state-feedback, 
with the remaining feedforward gain matrix G c derived from the system 


x(t) - K ff y re f(t) 


(FF GDcHs) GC c G 


x(t) 

x c (t) 
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or, in solving this (usually) nonsquare system for the steady-state response, 
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with as previously defined. 
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6.4 Concluding Remarks 

In output feedback design, consideration of closed-loop stability and disturbance rejection 
is done through careful definition of the input model. Under closed-loop transfer recovery, 
these design issues are first addressed under the setting of state feedback design. Many 
tradeoffs can be explored quickly using a state-feedback design, expecially when compared 
to the long turnaround associated with the direct optimization of an output-feedback con- 
troller. 

Direct optimization of an output-feedback controller may not include determination of 
command feedforward gains. Several useful approaches are given to compute feedforward 
gains for an output-feedback controller optimized using the CLTR approach. 
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Chapter 7 


Control Design for a 
High-Performance Rotorcraft via 
CLTR 


7.1 Introduction 

High-performance rotorcraft controller design is characterized by having to compensate 
for both the longitudinal and lateral dynamics in a single controller. Rather than having 
these modes largely decoupled in the natural state, one often needs to incorporate mode 
decoupling as a part of the controller design, especially when the system model includes 
high-frequency dynamics from the rotor flap and lag states. In addition to the usual design 
considerations of stability augmentation, a high-performance rotorcraft such as the UH-60 
would also require good bandwidth and decoupling in the heave, yaw rate, pitch, and roll 
command responses. 

The initial design step in CLTR is to formulate and synthesize a state-feedback controller 
that addresses all the design specifications. The synthesis consists of selecting appropriate 
design weights in the usual LQ performance index. For the UH-60 rotorcraft we construct 
an integral control, and incorporate a desired set of command input models. With integral 
control on all four important quantities — pitch, roll, yaw rate, and heave rate, we achieve 
the desired attitude and rate control of the respective output command variables as specified 
in ADS-33C for a rotorcraft. 

A series of output-feedback designs are then developed from this state-feedback de- 
sign using the procedure of Closed-Loop Transfer Recovery (CLTR). Initially, a Luenberger 
observer-based controller is designed by minimizing the difference between its closed loop 
response and the target closed-loop transfer function achieved with the state-feedback con- 
troller. The design minimization problem can be solved analytically through the solution of 
an algebraic Riccati equation. Note that in the recovery design process, we no longer need 
to iterate on the design performance index. With varying degrees of success in maintain- 
ing satisfactory closed-loop performance and robustness, other output-feedback designs are 
then derived from this controller using model reduction techniques. A numerical procedure 
is also used to design a CLTR controller following the problem description in Section 7. 1 1 . 
The controller is synthesized using one of these reduced-order designs as a starting point. 

On the other hand, one can also synthesize a reduced-order controller directly using the 
technique described in Chapter 2. Conveniently, one of the reduced-order designs devel- 
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oped from the Luenberger observer-based controller can be used as a starting point in the 
direct optimization. The resulting optimized design will be compared to the CLTR-based 
controllers. 

Design performance and robustness are evaluated based on the specifications defined in 
ADS-33C for a rotorcraft. 

7.2 The UH-60 Model 

The UH-60 rotorcraft model considered in this study is derived from a nonlinear simulation 
model linearized about two flight conditions: hover (nominal with 1-knot forward velocity) 
and a 15-knot forward velocity condition. In both cases the gross weight is set at 16,800 lbs, 
the rotor speed at 27 rad/s, and the air density at 0.002030 slug/cubic foot. The linearized 
open-loop model is of 31 st order. A listing of the rotorcraft states is given in Table 7.1. 


Table 7.1: UH-60 Rotorcraft States 


V, Q, r 

Body-axis attitude rates (deg/s) 

u , v> w 

Body-axis velocities (ft/s) 

Ox, 0y, e z 

Euler angles (deg) from terrestrial to 
rotorcraft axes (pitch, roll, yaw) 

x, y, Z 

Inertial positions (ft) 

$0, $D, 

PlS, PlC 

Rotor flapping rates (deg/s) 

Co, C D, 
Cis, Cic 

Rotor lead-lag rates (deg/s) 

Po, Pd, 
Pis, Pic, 

Rotor flapping angles (deg) 

Co, C D, 
Cis, Cic, 

Rotor lead-lag angles (deg) 

Ao, Ais, 
A ic 

Inflow velocities (1/sec) normalized 
by the rotor radius (26.83ft) 


The open-loop linear models in both the hover and 15-knot forward flight conditions 
are mildly unstable. The unstable poles represent a phugoid-like response in the front/side 
velocity coupled with the pitch and roll responses. In addition, there are modes that are at 
or near the origin corresponding either to a pure integration (for instance, the yaw variable 
derived from yaw rate) and to a lightly damped motion in the roll axis. Sensor outputs of 
the model consist of the body angular rates p, q, and r, the body velocities u, v, and w, the 
attitude angles roll (p, pitch 9, and yaw ip, and the inertial positions x, y, and z. 

In both flight conditions, the linear models have a pair of defective degenerate eigenvalues 
at the origin. In forward flight, the yaw angle couples into the lateral yi at displacement by 
way of forward velocity with y Jat (f) = V 0 ip(t ) included as part of the system model equations. 
(This is also true in the hover flight condition where there is a small, but not insignificant 
forward velocity) . These degeneracies within the system matrix and the tendency to produce 
overlapping modes within the CLTR design procedure are major concerns in the numerical 
solution of the design optimization. These concerns are addressed directly with the robust 
algorithm presented in Section 3.3. 
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The four control actuator inputs are <5 0 , S 3 , S c , and 6tr denoting the collective (main 
rotor), the main rotor sine, the main rotor cosine, and the tail rotor pitch inputs respectively. 
Using these controls, one should be able to achieve a decoupled set of ideal command 
responses in heave, roll, pitch, and yaw rate. 


7.3 Actuator and Sensor Delays 


Only the actuator delays and the sensor delays for the command variables 9, <j>, r and z 
(Figure 7.2) are present in the state-feedback design. Both the actuator and all the sensor 
delays are accounted for in the output-feedback designs. These were set uniformly to 60 ms. 
Because of the sensor delays, the system has a set of nonminimum-phase zeros and it is 
impossible to recover state- feedback design properties with an output-feedback design. 

This delay is modeled as a first-order Pade approximation model. In the Laplace domain, 
it is given by 

1 _ !* 

Tdelay{s) = , , 

1 + 2 


or, in the time domain, 


*(0 = + y in (t) 

4 

Voutit) = J,x{i) - y in (t) 

A higher-order Pade approximation model for the time-delay could have been used; however, 
this would lead to a prohibitive number of delay states due to the large number of sensor 
outputs (a total of 12) in this design problem. 

7.4 ADS-33C Requirements and the Ideal Response Model 

Performance of a controller design can be evaluated based on a subset of tests defined in 
the document ADS-33C [35]. This subset of requirements is applicable to rotocraft models 
linearized about specific flight conditions, and falls into the categories of bandwidth, attitude 
quickness, cross coupling, and gust response. 

Desired bandwidth properties of the closed-loop command responses are used to define 
low-order idealized response models for each of the command variables. The state-feedback 
design is derived from the minimization of the error between the actual rotorcraft responses 
and the ideal model responses. Such a procedure provides designers a systematic way of 
attaining adequate response in the respective channels without imposing excessively high 
bandwidth upon the others. Also, similar to CLTR, the need to adjust the criterion weights 
to optimize the design is not as severe as it would be with a direct design approach. The 
ideal model responses are associated with pitch, roll, yaw rate, and heave commands, and 
form the dynamics of the feedforward controller in the final design. 

Bandwidth requirements for both the forward and the hover flight conditions are similar, 
hence one feedforward command model is used for both conditions — with the exception that 
for turn coordination, a command for yaw rate has to be included with any roll command 
in forward flight. Otherwise, all commands are completely decoupled in these ideal models. 

Based on the ADS-33C bandwidth requirements in Figure 1(3.3) of [35], pitch, roll, 
and yaw rate responses are idealized as critically damped second-order responses, while 
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heave motion is patterned after a first-order response (refer to Table 4(3.3) in the ADS-33C 
document). A first-order approximation to a time delay of 0.1 s is included in each of the 4 
channels. As a result, the ideal model has a total of 11 states. The natural frequency of the 
pitch response is 3 rad/s, that of the roll response is 5 rad/s, and the yaw rate response has a 
bandwidth of 4 rad/s. Damping of the critically damped second-order responses is perturbed 
slightly away from 1 to avoid degenerate eigenvalues in these uncontrollable modes. These 
modes can cause numerical difficulties in the eigenvector decomposition method of the LQ 
synthesis. 

According to Figure 2(3.3) of ADS-33C, the command bandwidth is determined from the 
frequency response of the transfer function between the command input and the commanded 
output where the phase crosses 135 degrees. Targeted command bandwidths are set at 
4.33 rad/s for pitch response, 6.08 rad/s for roll response, and 5.27 rad/s for yaw rate 
response. Associated with each is a phase delay of 0.06s for pitch, 0.055s for roll, and 0.057s 
for yaw rate. While these are all Level I characteristics, the bandwidth is not excessively 
high. The heave response is set to a time constant of 3 s and a delay of 0.1s, which correspond 
to moderate Level I characteristics. These ideal responses are shown in Figure 7.6. 

Simplicity of the control synthesis based on ideal model matching allows designers to 
address the bandwidth problem as well as other design considerations. One area of con- 
cern is the attitude quickness. It is defined as the ratio between the peak rate and the 
angular change in an attitude step command, and is an evaluation criterion quantified in 
Figure 4(3.3) of the requirements [35]. These requirements are subject to interpretation 
when testing with a linear model. Because the criterion curve for ratio of maximum rate to 
attitude change is strongly dependent on the attitude command level, the most stringent 
case (at a 5° commanded change) is considered. To satisfy these criteria, ideal models with 
a pitch ratio of 1.1, and one with a roll ratio of 1.8 have been selected that correspond to 
average Level II responses in the target acquisition and tracking mode. To further improve 
these responses, one would have to increase the bandwidth further, and this would place 
a too stringent requirement upon any state or output-feedback designs in terms of control 
bandwidth. It is further anticipated that these ratios would be accentuated in a more re- 
alistic nonlinear simulation test, hence they are kept unchanged at the selected Level II 
characteristics. 

In the ADS-33C requirements described in items (3.3.9.2) and (3.4.4.2), pitch to roll cross 
coupling (and vice versa) is defined as the peak response in one attitude variable due to a step 
command in the other. Additional cross coupling tests in forward flight involve measuring 
the peak pitch response normalized by vertical acceleration in a step heave command (item 
3.4.4. 1.1), and the peak sideslip to a l°-step command in roll (item 3.4.6. 2). Additional 
hover condition tests examine the peak response and envelope of oscillatory response of yaw 
rate from a step command in heave (item 3. 3.9.1). 

Roll damping in forward flight is tested by measuring the envelope of the oscillatory roll 
response due to the roll command (item 3.4. 6.1). While neither pitch nor yaw damping are 
directly specified, the damping constants for the pitch, roll, and yaw rate responses are all 
set to near 1 in the ideal models. 

Additional hover condition requirements involve yaw and yaw rate responses to gusts 
(item 3. 3. 7.1). Tests for yaw responses to gusts are somewhat involved. In one test, one 
applies a steady 25-knot headwind, followed by a sudden 10-knot side gust. In the other 
test, the wind inputs are applied in the reverse order: a steady 25-knot side wind and then 
a sudden 10-knot head gust. 
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7.5 State Feedback Designs 


With the idealized output responses defined in the ideal models, the design procedure for 
the state-feedback designs follows somewhat the concept of loop transfer recovery — on the 
commanded responses. No reasonable ideal model appears to exist for the gust responses. 
By minimizing the difference in the commanded responses, the use of a pole attractor 
formulation for desired closed-loop stability is no longer needed. A further advantage is the 
anticipation factor introduced in the synthesis through the inclusion of the actuator delay 
in the design. A block diagram of the design set-up is given in Figure 7.1. 



Figure 7.1: Synthesis Model for Full-State Design with Integral Controls 


From our previous work [34], we recognize early on the need for integral control on 
the commanded error quantities. The purpose is to provide steady-state command track- 
ing to step changes and incorporate design robustness to constant disturbances. Putting 
penalty weights on the integral of the commanded errors ensure that the integral modes 
are detectable from the LQ design cost function. The penalty weighting matrices Q and R 
differ significantly between designs for the hover condition and forward flight. The control 
penalty is evenly divided among the 4 controls with R = I (identity matrix). The criterion 
variables consist of A 6 (pitch error), A <f> (roll error), Ar (yaw rate error), A z (heave error), 
and integral of the errors f AO, JAcfi, /Ar, and /Ai. In the hover condition, a criterion 
weight of 400 is used on all the variables with Q = 4007. In forward flight, however, a 
trade-off between performance and stability of an oscillatory mode around 18.5 rad/s is 
needed. This mode is excited primarily by the roll command. Sideslip control is of course a 
problem in forward flight, but an attempt to control it via a penalty on the sideways accel- 
eration v is not fruitful. Design effort involved in the selection of the penalty weights for the 
forward flight condition is considerably greater than that for the hover condition. The final 
weighting factors are given in Table 7.2. Robustness analysis in terms of single-loop gain 
and phase margins is also performed for both the control and sensor loops. Performance is 
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Table 7.2: Design We ights on the Criterion Variable s in the LQ Synthesis 


Criterion 

Variable 

Hover 

Forward 

Flight 

JA9 

400 

500 

I A0 

400 

500 

JAr 

400 

80 

JAz 

400 

100 

AQ 

400 

300 

Acp 

400 

400 

Ar 

400 

400 

Az 

400 

500 


satisfactory based on the ADS-33C requirements. Results are shown in Tables 7.7 to 7.10, 
and Figures 7.7, 7.8, 7.9, and 7.10. 

Since there is no penalty weighting on ip, x, y , or z, and that these states are associated 
with poles at the origin, they will not be detectable in the LQ synthesis. To avoid numerical 
difficulties in the solution of algebraic Riccati equation, these states have been removed 
from the state-feedback synthesis model. However, in the output-feedback designs, these 
additional sensor outputs would provide improved observability on the integrals of ip, u, 
v, and w. This information is especially helpful in steady-state command following for the 
variables r and z. 



Figure 7.2: Synthesis Model for Full-State Design with Sensor Delays Added 


With a state-feedback control structure as shown in Figure 7.1, all the state information 
is assumed available, uncorrupted by noises and not subject to sensor delay. However, for 
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proper control anticipation, sensor delays are taken into account in those sensor outputs that 
are used to form the integral errors. Hence, the state-feedback gains used in the subsequent 
CLTR observer-based designs are derived from the system shown in Figure 7.2. Performance 
of the resulting closed-loop state-feedback system defines the baseline results for comparison 
with subsequent output-feedback designs from loop transfer recovery or direct optimization. 
The integrators in the integral control of the command errors are clearly more a part of 
the feeback controller than of the plant model, and the state-feedback gains on the integral 
states can be retained in the reduced-order observer-based designs. 


7.6 Closed-Loop Transfer Recovery 

Closed-loop transfer recovery (CLTR) offers a useful framework for the synthesis of output- 
feedback controllers starting from a satisfactory state-feedback control laws. Analysis of 
closed-loop recoverability using the Special Coordinate Basis (SCB) transformation reveals 
that the system between the command/disturbance inputs and the measurement outputs 
has 23 internal states, 19 output states, 8 stable invariant zeros (at @33.33 rad/s), and 12 
infinite zeros of order 1. Note that the system used in this analysis includes the sensor 
delay in the outputs to the integral control. In this SCB analysis, we include independent 
disturbance inputs to all the system states, with the exception of the output delay states 
(for pitch, roll, yaw rate, and heave), and the command error integral states. Again note 
that the states ip, x, y, and 2 are not involved in the LQ synthesis and hence are excluded 
from the analysis. 

Exact or asymptotic recovery is in general not possible. The existence of invariant zeros 
in the right half-plane due to time delay in both the inputs and outputs poses a constraint on 
the loop recoverability. However since these invariant zeros are far removed from the origin 
and hence will contribute only slightly to the non-recoverable error. The major difficult for 
loop recovery lies in the fact that the system is not left-invertible due to the presence of 
multiple disturbances entering into the majority of the system states. 

It should also be noted that a full-order observer-based design generally achieves less 
overall recovery than a reduced-order (Luenberger) observer-based design since the latter 
design has the inherent benefit of feedthrough terms from the direct feedback of the mea- 
surement outputs. 


7.7 Luenberger Design 

From early experiences with similar rotorcraft designs, a Luenberger design tends to achieve 
better recovery of the closed-loop LQ performance. With direct observation of the integral 
states and the feedforward ideal model states, a Luenberger design is also more appropriate 
than a full-order observer-based controller. Furthermore, due to the large model size, one 
can make use of the sensor outputs to minimize the number of states that need to be included 
in the dynamics of the reduced-order observer. This is conveniently done by replacing parts 
of the system states with the noise-free measurement outputs. 

Design of a Luenberger observer-based controller under the CLTR procedure is described 
in Sections 5.5 and 5.6. In terms of an // 2 -norm minimization of the recovery error matrix, 
one usually needs to solve a singular optimal control problem. Generally, a reduction in the 
recovery error is accompanied by an increase in the magnitude of the observer gain matrix. 
Due to the presence of only infinite zeros of order 1, the recoverable portion can be exactly 
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recovered by a reduced-order Luenberger observer-based controller. A simple diagram of 
the CLTR design formulation can be seen in Figures 7.3 and 7.4. 



Figure 7.3: Block Diagram for CLTR with Output-Feedback Controller C(s) 



l_ J 

Figure 7.4: Synthesis Model for CLTR with Observer-Based Controllers 
The overall results are satisfactory, with some degradation in command bandwidth and 
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relatively minor changes in performance (Tables 7.13 and 7.14). Robustness also remains 
acceptable (Tables 7.7 to 7.12). Some noticeable differences in the command responses are 
seen in Figures 7.7 to 7.10. 

The overall controller design consists of a 35 t,l -order Luenberger observer combined with 
an 1 1'' 1 -order feedforward controller, and an integral controller of A th order. This compares 
to the 47 t/l -order rotocraft model, which includes a 31 st -order model for the rotocraft dy- 
namics, a 4 th -order model for the actuator delay, and a 12 i/l -order model for the sensor 
output delay. 

To describe more precisely the controller structure, let’s designate the ideal model (i.e., 
feedforward controller) system as {A^ai, B ideal ideal, AdeaJ, and the state matrices of 
the Luenberger design as 


Aic 

Bie 

Clc 

Gic 

Pic 


Luenberger observer system matrix 
Control input distribution matrix 
Output distribution matrix (state to control) 
Sensor to observer input distribution matrix 
Sensor to control distribution matrix 


Of particular interest are the Gi c and P/ c matrices. Note that the inputs to the controller 
include not only the usual sensor outputs y 3 from the plant, but the integral states X] 
and the feedforward controller states x/f as well. One can then partition the matrix Gi c 


into 


Gf c , G Gf c F where G\ c = G[ c f = 0. The same partitioning applies to the matrix 


p lc = p£ } P[ c , Pf/\ ■ These matrices show up in the overall controller model as follows, 
depending on whether it is described in 


• an observer-like form: 



'/A 

Gideal 

0 


Dideal 

A c = 

0 

Aideal 

0 

, Pc = 

Pideal 


0 

0 

Me 


0 

C c = 

pi 

r lc 

Pi7 c lc 

]. 

Dc = | 

> Pfe 


C/5 

0 


°] 


0 

0 

Plc 


(7.1) 


with inputs yemd , ys , and u, or in 
• an alternate non-observer form: 



[ / A —Cideal 0 


Dideal 

C /5 

A c = 

O 

O 

, Pc = 

Bideal 

0 


[ B lc P[ c B l/: P/; F A k + B, cC fc J 


0 

Gl + BtePg, 

C c = 

pL p L F c fc ], 

Dc = 

0 



(7.2) 

with no explicit input from the control u. 


Note that the term C /5 distributes the sensor ouputs <p, 6, r, and 2 to the appropriate 
integral states. So far, for the purpose of robustness analysis, the observer-based form has 
been used. Proper selection of one form over the other would become important when we 
proceed to the reduction of the controller order. 
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7.8 Model Reduction on the Feedforward Controller 


While the need to reduce the number of states in the 35 t/l -order Luenberger observer design 
is evident, one could also reduce the order of the feedforward controller as well. Usually one 
tends to perform a reduction of the entire controller at once without regard to the internal 
structures of the resulting linear time-invariant controller. However, to apply the method 
of balanced truncation, one needs to partition out the integral states of the commanded 
errors. Furthermore, the dynamics of the feedforward controller are only coupled in one 
direction to the feedback portion of the controller, it will therefore be practical to do this 
reduction as a separate process. 

Finally, one needs to examine the question as to whether one should reduce the com- 
ponents of the compensator in an observer-like form or in the alternate form. Putting the 
controller into the latter form will add a dependency to the feedforward controller states, 
as seen in the Bi c p£ f term within the A c matrix of equation 7.2. This term represents 
the distribution into a portion of the controller related to the Luenberger observer design, 
in the same place that G F C F would have been. Output from this distribution would add 
more degrees of freedom to the model reduction. Even though there are more inputs, the 
observer-like form is actually simpler and therefore preferred. 

Given a controller in observer-like form, the subsystem associated with the feedforward 
controller dynamics has inputs w c = {0 C , <f> c , ipc, z c }, and produces as outputs the ideal 
responses used in the feedforward commands wff — {@FF> 'pFFt ipFF> zff}, along with 
the feedforward control uff- Hence, the resulting feedforward controller used for model 
reduction is represented by the following state matrices 


AfF = AidedL , BfF = Bidecd 


Cff = 


Gideal 


jdFF 

Me 


Dff 


Bideal 

0 


Remember that the output uff is associated with the distribution matrix P t FF ■ In the 
model reduction, we seek the lowest model order that produces the closest input/output 
characteristics for the given command inputs. The preferred technique is based on the 
balanced truncation method. Looking at the Hankel singular values as listed in Table 7.3, 
one could identify that a 7 th -order feedforward controller (i.e., a reduction of 4 states) 
appears satisfactory. The same conclusion applies to both flight conditions. 

One problem with this reduction technique is that the static relation between yemd and 
wff in the resulting reduced-order feedforward controller is no longer unity. Usually curve 
fitting is used to further improve the results of the truncation procedure. In this case, we 
can make use of the additional degrees of freedom in the direct feedthrough term which is 
not modified in the balanced truncation procedure. The curve fit is done using numerical 
optimization. It is found that with white-noise inputs, the integral of the truncation error 
must be penalized in order for the steady-state values of the reduced-order model to match 
those of the original model. The error was reduced to less than a percent of its original 
value, and the resulting feedforward controller matches its full-order counterpart well. 

In the feedforward ideal model, only diagonal terms are considered in the command 
input distribution matrix. In other words, there is no crossfeed from a pitch command to 
the ideal roll response. For the command- to-control portion, a full feedthrough matrix is 
allowed. The direct command-to-control submatrix is contained in the D c matrix of 


84 



Table 7.3: Hankel Singular Values for Feedforward Controller Model 



Hover 

15knot Forward 

1 


13.4440 

23.1376 

2 


13.3388 

13.6907 

3 


12.7467 

11.8742 

4 


11.2347 

6.9295 

5 


2.5749 

4.0842 

6 


1.7669 

1.5641 

7 


1.5000 

0.9651 

8 


0.0654 

0.1047 

9 


0.0478 

0.0476 

10: 

0.0438 

0.0263 

11: 

0.0195 

0.0202 


the controller in observer-like form. Namely, 


D c 


Dpf? 1 



In the alternate form we have the same matrix in D c plus an additional nonzero submatrix 
Bi c D l pf^ 1 within the controller matrix B c , 


B c = 


Dideal 


Bideal 
d ryideal 

tilcUpic 


Cis 


0 

Gi + BicPg 


7.9 Balanced Truncation of Luenberger Observer 


The observer part of the compensator has not only sensor inputs y s , but also the control 
inputs u as well. These control inputs could be eliminated by separately closing the loop 
around the compensator, but this has two harmful side effects. Firstly, it makes impossible 
to identify and separate out in the compensator portions that correspond to the feedforward 
controller, the integral control, and the observer. Secondly, the open-loop compensator may 
possess some unstable poles, making the balanced realization impossible. Thus, it is essential 
to retain a separate input from the control u in what one would call an observer-like form 
of the compensator. 

The observer portion of the controller model has the following state model description, 


x 0 = A 0 x 0 + B 0 


Vs 

u 


u' = C 0 x 0 + D 0 u 


where 


A 0 — Ai c , B 0 


Die Bu-\ > C Q — Cic , D a — P^. 0 . 


Examining the Hankel singular values of such a system (Table 7.6), we see that the model 
can possibly be reduced to a 25 t/l -order controller without noticable change in performance. 

Design analysis results confirm the fact that this truncated controller model performs 
almost as well as the full-order observer-based controller design. In most areas of robustness, 
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command response, and ADS-33C design requirements the reduced order design shows little 
difference from the full-order design. 

Combined with the reduced-order feedforward controller (7 states) and the integral 
controller (4 states) , the overall controller is reduced to 36 t/l order. 

Closer examination of the Hankel singular values reveals that further reduction in the 
order of the observer part of the controller might be possible leading to a 14 t/l -order ob- 
server. Analysis of the controller with a 14 ift -order observer indicates that the design is not 
satisfactory, although the resulting closed-loop system remains stable. In the next section, 
this reduced-order observer will be redesigned using numerical optimization for improved 
performance and robustness. 


7.10 Reduction of Controller Order using Numerical Opti- 
mization 


Note that the feedforward controller order has been reduced separately, and the integral 
gains are left unchanged, the remaining task in controller order reduction involves only the 
observer portion of the controller. As such, only this portion of the controller is formulated 
for optimization using the design tool SANDY. 

The controller state matrix A c matrix is allowed to have a tridiagonal structure which 
was derived from model reduction of the Luenberger observer. All of the elements in the 
matrix A c identified by the asterisk are design variables in the optimization. 

'** 000000000000 ' 
***00000000000 
0***0000000000 
00***000000000 
000***00000000 
0000***0000000 
_ 00000***000000 
Ac ~ 000000***00000 
0000000***0000 
00000000***000 
000000000***00 
0000000000***0 
00000000000 *** 
. 000000000000 **. 


Furthermore, elements of the remaining controller state matrices B c , C c> and D c are also 
choosen as design parameters. 


B c 

C c 

D c 


(12 sensor inputs) (4 control inputs) 


(4 control outputs) C[ c 


(4 control outputs) 



(4 control inputs) 


While the \A th -order observer-based controller arrived at by balanced truncation is inad- 
equate, it can however be used as a starting point from which one can recover the full-order 
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controller responses through numerical optimization. We define the design objective func- 
tion in the design algorithm SANDY to be the H 2 - norm of the error in transfer function 
between the original Luenberger observer and the above specified 14 t/l -order design. As in 
the optimization of the feedforward controller, it would be desireable to incorporate more 
design degrees of freedom into the lower-order observer structure without increasing the 
overall controller order. However, because the optimization already involves 368 variables, 
introducing additional degrees of freedom will have to be done carefully. One could for 
example allow a direct feedthrough term from the control input u to itself, but it is found 
to be unnecessary for a well-posed system. One could also have introduced a direct link 
from the dynamics of the feedforward controller and the integral control to the observer 
dynamics. However, this path has been partially fulfilled by the control input into the 
observer and the formulation would be redundant. 

In consideration of the alternate formulation to the observer-like form, it is found that 
closing the control loop about the compensator is not beneficial. While it could be argued 
that eliminating the explicit control input u will reduce the parameter count, a closer 
examination indicates that there is no significant advantage. 

Simultaneous optimization of parameters in both the feedforward controller and the 
Luenberger observer controller is more costly in terms of memory and computational time. 
Optimization of just the observer-based controller alone involves 368 parameters. Because 
the feedforward controller matches its full-order counterpart well, and because the results 
from separate optimization of the observer-based controller are satisfactory, the procedure 
involving simultaneous optimization is therefore not needed. 

The resulting design optimization improves significantly the controller performance. 
While the performance is close to the other output-feedback controllers (i.e., small rela- 
tive to the difference between the other output-feedback controllers and the state feedback 
controller), its robustness while still adequate is somewhat degraded. 


7.11 Numerical CLTR 

In general, one can also approach the closed-loop transfer recovery via direct numerical 
optimization. It should be emphasized that numerical CLTR is different than the procedure 
described in Section 7.10 related to controller order reduction. Here, one can theoretically 
start from any preferred size of the controller — the order selected is an arbitrary choice. 
The design objective is to determine the controller parameters by matching the closed- 
loop responses of the resulting output-feedback controller to those achieved under the state 
feedback for a given set of disturbance/command inputs. 

In the design optimization, the 25^-order controller (derived from the 14 th -order ob- 
server given in the previous section) provides a convenient starting place in this case due to 
its size, and its potential for further improvement. 

Initially, the input excitations are bandwidth-limited to 25 rad/s and fed into both the 
command and gust input channels. The resulting controller design exhibits poor robustness 
in the control paths. Additional noises bandwith-limited to 150 rad/s are then introduced 
into the control actuators (the high-bandwidth value is choosen to be above all the system 
modes). For the hover case, the intensity of the actuator noises is set at 0.666 compared 
to a unit intensity in all other input excitations. However, the intensity is increased to 1.0 
for the forward flight case (The results are found to be relatively insensitive to the choosen 
value). 
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Included in the performance index are the errors in the command responses in roll, 
pitch, yaw rate, heave, and their integrals. For robustness, it is necessary to also include 
the difference in control responses in the performance index. Several design iterations are 


Table 7.4: Objective Weights in Numerical CLTR 


Criterion 

Hover 

Forward 

Variable 


Flight 

Roll Residual A <j) 

400 

300 

Pitch Residual Ad 

400 

400 

Yaw Rate Residual Ar 

SI 

400 

Heave Residual A i 

800 

500 

/ A0 

700 

500 

JAd 

400 

500 

f Ar 

400 

80 

JAz 

100 

100 

Collective Residual A<5o 

1 

1 

Main Rotor Sine Residual AS s 

1 

1 

Main Rotor Cosine Residual A6 C 

1 

1 

Tail Rotor Residual AStr 

1 

1 


usually needed in the set up and selection of the weightings in this strongly non-recoverable 
case. Table 7.4 lists the weightings used— they are similar to, and based on the state- 
feedback performance objective. The optimization results show that no further improvement 
can be made to the design obtained in Section 7.10. 

7.12 Direct Optimization 

In this section ,we examine design results obtained from the direct optimization of an 
objective function instead of applying the CLTR procedure. Here, the controller is designed 
by minimizing the difference between the ideal command responses and the actual responses 
along with terms involving control penalty. Initial values of the criterion weights are taken 
from the respective Q and R matrices defined in the LQ synthesis. The controller structure 
is chosen to be the same as that in the last section to further explore the potential for 
optimization beyond the current designs having the same order and structure. The controller 
order is found to be reasonably small for design implementation. The numerical CLTR 
results of Section 7.10 are used to provide initial design guess in the direct optimization. 

As with the numerical CLTR, the system is excited by unit amplitude command and gust 
inputs as well as by fictitious actuator noises. The command and gust inputs are modelled 
as outputs of first-order filters with a bandwidth of 25 rad/s excited by white noises of 
unit intensity. The actuator noises are derived from first-order filters with a bandwidth of 
150 rad/s. While the actuator noises are kept at unit amplitude in the forward flight case, 
they are adjusted somewhat in the hover case: 0.2 for the collective control, and 0.5666 for 
the other control actuators. 

In the hover case, considerable improvement can be achieved. The weighting matrices 
are initially set to the values defined in the LQ synthesis, and undergo some adjustment 
through several design iterations. The final values can be seen in Table 7.5. A slight change 
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Figure 7.5: Controller Design Structure in Direct Optimization 


Table 7.5: Objective Weights in Direct Optimization 


Design 

Hover 

Forward 

Variable 


Flight 

Roll Residual A 

500 

300 

Pitch Residual AO 

500 

400 

Yaw Rate Residual Ar 

400 

400 

Heave Residual A i 

800 

500 


750 

500 

I AO 

650 

500 

f Ar 

400 

80 

JAz 

300 

100 

Collective Control So 

1 

1 

Main Rotor Sine Control 6 S 

1 

1 

Main Rotor Cosine Control 8 C 

1 

1 

Tail Rotor Control Str 

1 

1 “ 
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in the forward flight design is found, but the difference is so slight that there is no noticeable 
changes in the design results. It would be difficult to characterize the difference in results 
for the forward flight, except that direct optimization may offer a slight edge in the overall 
robustness. Most of the improvement in the hover design is in the roll bandwidth (see 
Table 7.7). The improvement is gained at the expense of sensor robustness (Table 7.13). 

7.13 Conclusions 

Applicability of the CLTR technique has been demonstrated on the control of the UH-60 
rotorcraft for two flight conditions. The results are presented for the hover and 15-knot 
forward flight conditions. They are found to be satisfactory in terms of ADS-33C require- 
ments. It turns out that the feedforward controller dynamics at the two flight conditions are 
nearly the same, but a low-order set of ideal response dynamics should be feasible for any 
flight condition. With the given ideal model dynamics, a reasonable state-feedback design 
can be synthesized. 

Not all the requirements in ADS-33C can be tested in the regime of linear models. 
Generally the controllers are found to meet the requirements for Level I qualities, except, as 
noted, for pitch and roll rate responses. The 14 th -order Luenberger observer-based controller 
for hover does not meet Level I characteristics for roll bandwidth, and the result has been 
improved via direct optimization. This is not indicative of the limitations in the CLTR 
design procedure. It is believed that an improved state-feedback design for the hover case 
would have led to better output-feedback controllers via CLTR. 

The proposed design technique based on CLTR should be applicable to other flight 
conditions. Evaluation of the CLTR design process for the hover and the forward (15 knot) 
flight conditions shows that most of the design tradeoffs occur in the synthesis of the state- 
feedback design. Once the problems of stability, performance and robustness have been 
solved in the LQ synthesis, recovery of these closed-loop properties becomes a much more 
routine process under CLTR. Although the procedure of numerical CLTR requires the 
tailoring of objective weighting matrices for an optimum design at each flight condition, an 
accurate selection of the weighting matrices is not significant to the overall results. 

The order of the overall synthesis model in a CLTR design using numerical optimization 
is nearly prohibitive. Some shortcuts are performed to reduce the size of the overall CLTR 
system. A reduced-order (7 state) feedforward controller is used to drive both the state 
feedback (i.e., target) closed-loop model and the nominal open-loop model. The state- 
feedback system with 4 sensor output delay states is 3b th order, while the open-loop plant 
with a full set of sensor output delays is 47 th order. With integral control (4 states) for both 
state- and output-feedback, and a 14 t/l -order observer, the overall dynamics quickly add up 
to a synthesis model of 111 states. An optimization on such a system takes substantial 
computing time, and is prone to have defective eigenvalues throughout the optimization 
run. The design results rely heavily on the reliable algorithm developed in Chapter 3. 

The procedure in direct optimization design is facilitated considerably by the results from 
CLTR. The initial choice of objective weights can be developed from the LQ synthesis, and 
the initial controller can be taken from the CLTR design. In fact, the choice of controller 
order would not have been possible without the CLTR design followed by a model reduction 
procedure. One can conclude that design based on direct optimization would have entailed a 
much more tedious process if not for the useful insights developed from the CLTR procedure. 
The CLTR process becomes truly an integral part of the direct optimization design. 





Table 7.6: Hankel Singular Values, Luenberger Observer 



Hover 

15knot Forward 

1: 


19.6542 

22.3400 

2: 


16.1352 

18.2805 

3: 


14.3219 

15.3796 

4: 


11.1515 

10.5320 

5: 


9.9974 

9.0218 

6: 


6.5465 

7.0539 

7: 


4.7171 

6.7871 

8: 


3.4355 

4.8661 

9: 


3.2165 

3.3467 

10: 

2.7132 

3.2328 

11 


2.1561 

2.4144 

12 


1.8829 

1.8143 

13 


1.6118 

1.2724 

14 


1.1177 

0.7836 

15 


0.5927 

0.6148 

16 


0.4878 

0.4694 

17 


0.3356 

0.3540 

18 


0.2910 

0.2724 

19 


0.2061 

0.1688 

20 


0.1795 

0.1533 

21 


0.1721 

0.1211 

22 


0.1270 

0.1111 

23 


0.1071 

0.0884 

24 


0.0967 

0.0551 

25 


0.0190 

0.0495 

26 


0.0003 

0.0030 

27 


0.0002 

0.0025 

28 


0.0000 

0.0005 

29 


0.0000 

0.0000 

30 


0.0000 

0.0000 

31 


0.0000 

0.0000 

32 


0.0000 

0.0000 

33 


0.0000 

0.0000 

34 


0.0000 

0.0000 

35 


0.0000 

0.0000 


Table 7.7: Single-Loop Sensor Robustness Properties, Hover 


Design 

Roll 

Gain 

Phase 

Pitch 

Gain 

Phase 

Heave 

Gain Phase 

Yaw Rate 
Gain Phase 

Full State 
-j-Integrals 

6.13db, 

<-40db 

52.36 

8.84db, 

<-40db 

48.85 

+7.15db, 

<-40db 

58.48 

+5.97db, 

<-40db 

57.39 

Luenberger 

Observer 

+5.64db, 

<-40db 

49.92 

+6.88db, 

<-40db 

49.75 

+5.01db, 

<-40db 

54.80 

+6.31db, 

-16.04db 

44.75 

25th Order 
Luenberger 

+6.43db, 

<-40db 

49.93 

+6.60db, 

<-40db 

49.81 

+5.39db, 

<-40db 

50.17 

+6.73db, 

<-40db 

44.77 

14th Order 
Luenberger 

+3.98db, 

-12.41db 

45.00 

+5.60db, 

<-40db 

46.68 

+5.44db, 

<-40db 

52.98 

+5.24db, 

-13.78db 

43.27 

14th Order 
Direct Opt. 

+4.03db, 

-6.89db 

34.95 

+5.98db, 

<-40db 

43.71 

+6.29db, 

-28.96db 

58.93 

+5.07db, 

-29.20db 

41.95 
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Table 7.8: Single-Loop Actuator Robustness Properties, Hover 


Design 

Collective 
Gain Phase 

Sine 

Gain 

Phase 

Cosine 

Gain Phase 

Tail Rotor 
Gain Phase 

Full State 
-l-Integrals 

>40db, 

<-40db 

65.00 

>40db, 

<-40db 

80.63 

>40db, 

<-40db 

70.80 

>-J-40db, 

<-40db 

91.55 

Luenberger 

Observer 

A V 

1 hti. 

2 CL 

65.00 

>40db, 

<-40db 

84.59 

>40db, 

<-40db 

73.97 

>+40db, 

<-40db 

91.55 

25th Order 
Luenberger 

>40db, 

<-40db 

65.00 

>40db, 

<-40db 

80.63 

>40db, 

<-40db 

70.79 

>+40db, 

<-40db 

91.55 

14th Order 
Luenberger 

>40db, 

<-40db 

48.40 

>40db, 

<-40db 

74.97 

>40db, 

<-40db 

59.89 

>+40db, 

<-40db 

97.21 

14th Order 
Direct Opt. 

>40db, 

<-40db 

61.28 

>40db, 

<-40db 

57.46 

>40db, 

<-40db 

67.38 

>+40db, 

<-40db 

103.82 


Table 7.9: Multiloop Actuator Robustness Properties, Hover 


Design 

Independent 
+G.M. -G.M. 

P.M. 

<Z(I + G(s)) 

*(I + C-'(s)) 

Full State 

-foe db, 

-lOdb 

60° 

0.992 

0.693 

Luenberger Observer 

-foe db, 

-lOdb 

60° 

0.992 

0.693 

25th Order Luenberger 

-foe db, 

-lOdb 

60° 

0.992 

0.693 

14th Order Luenberger 

+12, 

-7 

40 

0.715 

0.571 

14th Order Direct Opt 

+12, 

-7 

40 

0.700 

0.607 


Table 7.10: Single-Loop Sensor Robustness Properties, Forward Flight 


Design 

Roll 

Gain 

Phase 

Pitch 

Gain 

Phase 

Heave 

Gain Phase 

Yaw Rate 
Gain Phase 

Full State 
-j-Integrals 

6.86db, 

<-40db 

49.89 

8.48db, 

<-40db 

47.70 

+6.16db, 

<-40db 

61.03 

+7.24db, 

<-40db 

61.20 

Luenberger 

Observer 

+5.49db, 

<-40db 

49.70 

+6.82db, 

<-40db 

44.99 

+6.89db, 

<-40db 

61.03 

+4.39db, 

<-40db 

44.94 

25th Order 
Luenberger 

+5.38db, 

<-40db 

49.36 

-f6.84db, 

<-40db 

44.99 

+6.89db, 

<-40db 

64.80 

+5.29db, 

<-40db 

43.91 

14th Order 
Luenberger 

+5.41db, 

<-40db 

44.70 

+6.56db, 

-28.61db 

43.88 

+6.48db, 

<-40db 

54.95 

+5.43db, 

<-40db 

42.01 

14th Order 
Direct Opt. 

+5.41db, 

<-40db 

44.70 

+6.85db, 

-28.61db 

43.88 

+6.48db, 

<-40db 

54.95 

+4.61db, 

<-40db 

42.01 
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Table 7.11: Single Loop Actuator Robustness Properties, Forward Flight 


Design 

Collective 
Gain Phase 

Sine 

Gain 

Phase 

Cosine 

Gain Phase 

Tail Rotor 
Gain Phase 

Full State 
-{-Integrals 


64.29 

>40db, 

<-40db 

74.95 

>40db, 

<-40db 

74.90 

>+40db, 

<-40db 

90.34 

Luenberger 

Observer 


64.95 

>40db, 

<-40db 

74.77 

>40db, 

<-40db 

75.00 

>+40db, 

<-40db 

90.34 

25th Order 
Luenberger 


64.28 


74.84 

>40db, 

<-40db 

75.00 

>+40db, 

<-40db 

90.34 

14th Order 
Luenberger 

>40db, 

<-40db 

45.93 

>40db, 

<-40db 

54.34 

>40db, 

<-40db 

57.23 

>+40db, 

<-40db 

90.44 

14th Order 
Direct Opt. 

>40db, 

<-40db 

45.93 

>40db, 

<-40db 

54.94 

>40db, 

<-40db 

59.51 

>+40db, 

<-40db 

90.44 


Table 7.12: Multiloop Actuator Robustness Properties, Forward Flight 


Design 

Independent 
+G.M. -G.M. 

P.M. 

<L(I + G{s)) 

<t(I + G-\s)) 

Full State 



■sum 

1.000 

0.632 

Luenberger Observer 



K&H 

1.000 

0.632 

25th Order Luenberger 






14th Order Luenberger 

+12, 

-9 




14th Order Direct Opt 

+12, 

-9 

40 

0.709 
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Table 7.13: Design Evaluation, Hover 


Test Result 

ADS-33C 

Target 

Full State 
Feedback 

Luenberger 

Observer 

25th Order 
Luenberger 

14th Order 
Luenberger 

14th Order 
Direct Opt. 

Pitch 

Bandwidth 

Fig la,c 

3.344 

2.642 

2.648 

2.699 

2.510 


Phase Delay 

(3.3) 

0.147 

0.187 

0.186 

0.184 

0.182 

Roll 

Bandwidth 

Fig lb,d 

5.333 

3.626 

3.649 

2.819 

3.485 


Phase Delay 

(3.3) 

0.080 

0.120 

0.120 

0.086 

0.104 

Yaw 

Bandwidth 


5.311 

KI1 H 

4.116 

4.240 

4.130 

Rate 

Phase Delay 

Hs&jfli 

0.061 


0.105 

0.102 

0.113 

Heave 

Time Constant 

< 5.0 

2.981 


2.856 

2.816 

2.876 


Delay (s) 

< 0.2 

0.114 

0.169 

0.173 

0.171 

0.195 

Pitch 

Rate from 1° step 


1.160 

1.177 

1.175 

1.234 

1.220 


from Roll 

wbsm 

0.012 



0.016 

0.007 

Roll 

Rate from 1° step 

> 2.4 

1.813 



1.914 

2.020 


from Pitch 

< 0.25 

0.009 

0.011 


0.022 



from side gust 


0.264 

0.292 

0.292 

0.297 


Yaw 

from head gust 


0.264 

0.292 

0.292 

0.297 

0.283 

Rate 

peak, from heave 


5.8e-4 

8.2e-4 

9.1e-4 

65. le-4 

59.7e-4 


oscil., from heave 


5.8e-4 

8.1e-4 

11.8e-4 

13.2e-4 

45.7e-4 


for 1 ° Yaw change 

| > 2.4 | 

2.800 

2.773 

2.835 

2.861 

3.430 


Table 7.14: Design Evaluation, Forward Flight 



Test 

ADS-33C 

Full State 

Luenberger 

25th Order 

14th Order 

14th Order 


Result 

Target 

Feedback 

Observer 

Luenberger 

Luenberger 

Direct Opt, 

Pitch 

Bandwidth 

Fig 1 

3.741 

3.301 

3.296 

3.325 

3.324 


Phase Delay 

(3.4) 

0.150 

0.190 

0.189 

0.180 

0.180 

Roll 

Bandwidth 

BESS 

5.779 

4.825 

4.810 

4.844 

4.844 


Phase Delay 

mmM 

0.176 

0.136 

0.136 

0.115 

0.115 

Yaw 

Bandwidth 

Fig 8 

5.329 

4.059 

4.063 

4.244 

4.244 

Rate 

Phase Delay 

(3.4) 

0.061 

0.105 

0.107 

0.112 

0.112 

Heave 

Time Const. 

<5.0 

2.957 




2.836 


Delay (s) 

<0.2 

0.115 




0.178 

Pitch 

from vertical accel. 


0.00779 

0.00779 

0.00867 

0.00822 

0.00822 


from Roll Cmd. 


0.002 

0.002 

0.002 

0.006 

0.006 


Rate from 1° step 

> 2.4 

1.860 

1.847 

1.839 

1.846 

1.846 

Roll 

Osc. from 1° step 

Fig 5(3.4) 

0.00223 

0.00217 

0.00151 

0.00172 

0.00171 


from Pitch Cmd. 

<0.25 

0.014 

0.013 

0.013 

0.012 

0.012 

Sideslip 

from 1° Roll step 

Fig 6(3.4) 

0.110 

0.110 

0.111 

0.105 

0.105 
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State Feedback 

Luenberger Observer 

25th Order Luenberger 

14th Order Luenberger w/SANDY 

— 14th Order Direct Optimization 
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0.02 



State Feedback 

Luenberger Observer 

25th Order Luenberger 

14th Order Luenberger w/SANDY 
— 14th Order Direct Optimization 



97 



State Feedback 

Luenberger Observer 

25th Order Luenberger 

14th Order Luenberger w/SANDY 
— 14th Order Direct Optimization 
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State Feedback 

Luenberger Observer 

25th Order Luenberger 

14th Order Luenberger w/SANDY 

14th Order Direct Optimization 
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Chapter 8 

Future Work 


8.1 Evaluation of Worst-Case Method for Loop Transfer 
Recovery 

It is perceived that minimizing the worst-case difference between the closed-loop responses 
of a state feedback design and a those corresponding to a low order output feedback con- 
troller case to a given range of disturbances would produce a more acceptable controller 
than from the H 2 norm of the error. Given reasonable recovery, such a controller would be 
less prone to extreme variation from the state feedback characteristics. This would include 
loop transfer properties used in robustness measures. There exist several good test problems 
where this design algorithm can be exercised. 

8.2 Completion of the Hybrid Algorithm 

The decomposition of the system matrix into a partition of non-defective eigenvalues and a 
partition of defective eigenvalue blocks, as detailed in Appendix A, needs further evaluation. 
A thorough test must preceed the conversion of the algorithms for computing X and M to 
the hybrid form. 

The hybrid procedures for evaluating X and M., given in Appendix B, are more intri- 
cate than either the robust or diagonal forms due to the presence of cross terms. These 
cross terms are integrals containing both non-defective eigenvalues and defective eigenvalue 
blocks. The additional numerical complexity was relatively slight for the X calculation 
where simple closed form solutions could be derived to reduce the complexity of the cross 
terms. Explicit handling of every cross term in a custom formula may yield rewards in 
speed and memory usage. However, in the case of M, the problem is much more difficult. 
For example, for the double integral containing the defective eigenvalue block Wc 

r e x ' {v - s) B l i 2 e Wcv V2ije x ’ s dsdv , 

Jo Jo 

it may still be as efficient to simply apply the robust algorithm using these arguments. 

8.3 Eigenvalue/Damping Constraints 

A tidy design is often achieved in pole placement. While exact pole placement is not 
possible within the context of numerical optimization, building a cost function based on 

>age—ML . INTENTIONALLY BLANK 
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pole attractors [3] or forcing designated poles into a selected region [6] would seem to 
be. Creating a cost function that includes eigenvalue constraints and yet is robust to the 
degenerate mode condition is a subject for further reseach. 


8.3.1 Current Method 

The current method for eigenvalue constraint depends on the eigenvalue — eigenvector de- 
composition, and thus will not work when degenerate eigenvalues exist in the closed-loop 
system model. 

The eigenvalue constraints and their gradients with respect to the controller design 
parameters are defined as follows. Given a system matrix A we let 


Avi — AjUi (8.1) 

where A j and v t are the eigenvalues and eigenvectors, respectively, of the system matrix A. 
We define a simple constraint function that is zero when the eigenvalues are in the desired 
stability region, or a nonzero value equal to the distance away from the desired stability 
boundary. One such function is 


Ja 


1 

— ) ' {max(ai ffmon 0 )} 
1 i=l 


This would reflect a cost function on real (or real part of) eigenvalues greater than a certain 
(T-max- It would be zero if all eigenvalues were less than a max- 

In addition to constraining the real part, one usually needs to constrain the damping 
ratio 


<i = - 



to be larger than a certain value (say Cmin)- Consider 


1 71 

Jc = » E {max(<7i sing + [oleosa, 0)} 2 

^ i= 1 


where cos a = ( min and sin a = Jl - Cmtn- Each summand of this function will be zero 
if the corresponding eigenvalue satisfies the constraint; otherwise, it will be equal to the 
square of the distance from the offending eigenvalue to the damping line. 

The above constraint functions are differentiable and their gradients are useful for numer- 
ical optimization, especially for the nonlinear programming algorithm NPSOL [21]. Given 

1 n 

J<x = — ) ] {l7l(lx((Ti (Tmaxi 0)} 


we have 



Also, the damping constraint gradient is 


dJ^ 

dp 


E 

i=l 


max{(Ji sin a + |wi| cos 


n\ ( d<Ji 

0) (a? 


sine* + sgn(ui )%^ cos 

dp 


“)} 
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For non-degenerate eigenvalues, these gradients are well behaved. 

Note that the constraint gradients are expressed in terms of the gradients of the closed- 
loop eigenvalues. Taking the derivative of equation (8.1) with respect to a parameter p, we 

~ ‘ . dv t dX i , x dVi 


dA 

—-Vi + A- 
dp op 


-Vi + A i- 

dp op 


or: 


(A-Xil) 


dvi d \ i dA 

-3— - -K~Vi = V, 

dp op op 


The objective is to find an equation without a ^ term. One can diagonalize: A = 
TAT -1 , where A is a diagonal matrix containing the eigenvalues of A. With this transfor- 
mation applied to the above equation, we have: 


Ai — A» 0 . . . . 0 

0 

0 Ai_i - Ai 0 

0 0 0 

. 0 Ai-)-i — Ai 0 

0 

0 . . . . 0 X n — Xi 


1 dvi 1 d X i 

T ir~T 

dp dp 


---I 


dA 
dp 1 


The i th equation from the above set shows a simple relation between ^ and without a 
term thus one can write: 


d Xj 

dp 


(r-'SH 


The notation (.)i denotes the i th row of the enclosed term. For complex eigenvalues, o + iu, 
the real and imaginary parts of the above are §| and respectively. Note that when there 
is a near degeneracy and the eigenvalues are nearly defective, T -1 will be badly conditioned. 


8.3.2 New Method of Eigenvalue Constraints 

There are two equally important aspects of these constraints. First is the mechanization of 
the constraints themselves. The second is the ability to identify those parts of the system 
matrix where the constraints should apply. One often has disturbable yet uncontrollable 
modes where the eigenvalues do not obey the constraints and would upset the optimization 
process if so included. Moreover integral poles are often formulated to serve a control or 
estimation purpose. 

It could be possible with the proposed method of decomposing the system matrix into 
eigenvalue blocks one could derive a means of calculating the eigenvalue derivatives for those 
blocks. All attempts so far have seemed preliminary. 

Finally it may become that the only really promising method for forcing designated 
eigenvalues to a given region will be the plant transformation methods of Kawasaki and 
Shimemura [6] or Bernstein and Haddad [8]. 


103 


104 


Bibliography 


[1] Stein, G. and Athans, M., “The LQG / LTR Procedure for Multivariable Feedback 
Control Design,” IEEE Transactions on Automatic Control , AC-32, pp. 105-114, 1987. 

[2] Brasch, F.M., and Pearson, J.B., “Pole Assignment Using Dynamic Compensation”, 
IEEE Transactions on Automatic Control, Vol AC-15, No.l, 1970, pp. 32-43. 

[3] Gordon, V.C., and Collins, D.J., “Multi-Input Multi-Output Automatic Design Syn- 
thesis for Performance and Robustness”, Journal of Guidance and Control, May- June 
1986, Vol. 9, No. 3, pp. 281-287. 

[4] Doyle, J.C., Glover, K. , Khargonekar, P.P. and Francis, B.A., “State-Space Solutions to 
Standard H 2 and // 00 -Control Problems,” IEEE Transactions on Automatic Control, 
Vol. AC-34, No. 8, pp. 831-847, 1989. 

[5] Doyle, J.C. and Glover, K., “State-Space Formulae for All Stabilizing Controllers that 
Satisfy an EToo-Norm Bound and Relations to Risk Sensitivity,” Systems & Control 
Letters, 11, pp. 167-172, 1988. 

[6] Kawasaki, N. and Shimemura, E., “A Method of Deciding Weighting Matrices in an 
LQ-Problem to Locate All Poles in the Specified Region”, IFAC Control Science and 
Technology (8th Annual World Congress) Kyoto, Japan, 1981. pp 481-486. 

[7] Kimura, H., “Pole Assignment by Gain Output Feedback”, IEEE Transactions on 
Automatic Control, Vol. AC-20, No.4, 1975, pp. 509-515. 

[8] Bernstein, D.S., and Haddad, W.M, “Controller Design With Regional Pole Con- 
straints”, IEEE Transactions on Automatic Control \ ol 37, No.l, January 1992, pp 
54-69. 

[9] Mills-Curran, W.C., “Calculation of Eigenvector Derivatives for Structures with Re- 
peated Eigenvalues”, AIA A Journal, Vol 27, No.7, July 1988, pp. 867-871. 

[10] Sannuti, P. and Saberi, A., “A Special Coordinate Basis of Multivariable Linear 
Systems-Finite and Infinite Zero Structure, Squaring Down, and Decoupling”, Int. 
J. Contr., Vol. 45, No. 5, pp. 1655-1704, 1987. 

[11] Chen, B.M., Saberi, A. and Ly, U., “Closed-loop transfer recovery with observer based 
controllers — Part 1: Analysis,” Proceeedings of the 1991 AIAA GNC Conference. Also, 
to appear in Control and Dynamic Systems: Advances in Theory and Applications. 

[12] Chen, B.M., Saberi, A. and Ly, U., “Closed-loop transfer recovery with observer based 
controllers — Part 1: Design,” Proceeedings of the 1991 AIAA GNC Conference. Also, 
to appear in Control and Dynamic Systems: Advances in Theory and Applications. 

[f. fifi£j^_WTENTiOMALLY BLAL’L 105 

PAKWHN6 PAGE BLANK NOT FILMED 



[13] Stoorvogel, A. A., ‘The Singular Hoo -Control Problem with Dynamic Measurement 
Feedback,” SIAM J. Control and Optimization , Vol. 29, No. 1, pp. 160-184, 1991. 

[14] Stoorvogel, A. A. and Trentelman, H.L., ‘The Quadratic Matrix Inequality in Singular 
//(» -Control with State Feedback,” SIAM J. Control and Optimization , Vol. 28, No. 5, 
pp.1190-1208, September 1990. 

[15] Levine, W.S. and Athans, M.,”On the Determination of the Optimal Constant Output 
Feedback Gains for Linear Multivariable Systems,” IEEE Transactions on Automatic 
Control, AC-15, pp.44-48, 1970. 

[16] Anderson, B.D.O. and Moore, J.B., Linear Optimal Control, Englewood Cliffs, Prentice 
Hall, Inc., New Jersey, 1971. 

[17] Makila, P.M., “Computational Methods for Parametric LQ Problems A Survey,” IEEE 
Transactions on Automatic Control, AC-32, No.8, pp.658-671, August 1987. 

[18] Ly, U., A Design Algorithm for Robust Low-Order Controllers, PhD. Thesis, Depart- 
ment of Aeronautics and Astronautics, Stanford University, November, 1982. 

[19] ElGhaoui, L., Carrier, A. and Bryson, A.E., “A New Performance Criterion for Worst 
Case Analysis and Design”, Journal of Guidance and Control, July- August, 1992, pp. 
953-961. 

[20] Mills, R.A., Robust Controller and Estimator Design Using Minimax Methods, PhD. 
Thesis, Department of Aeronautics and Astronautics, Stanford University, March 1992. 

[21] Gill, P.E., Murray, W., Saunders, M.A. and Wright, M.H., User’s Guide for NPSOL 
(Version 4-0): A Fortran Package for Nonlinear Programming. Technical Report SOL- 
86-2, Stanford university, January 1986. 

[22] VanLoan, C.F., “Nineteen Dubious Ways to Compute the Exponential of a Matrix”, 
SIAM Review 20, pp.801-836, 1978. 

[23] Golub, G.H. and VanLoan, C.F., Matrix Computations. Johns Hopkins University 
Press, pp 396-400, 1985. 

[24] Bender, C.M. and Orzag, S.A. Advanced Mathematical Methods for Scientists and En- 
gineers. McGraw-Hill 1978. pp 383-389. 

[25] Osborne, E.E. , ”On Preconditioning of Matrices”, Journal of the Association for 
Computing Machinery, Los Angeles, CA. pp 338-345. March, 1960. 

[26] Parlett, B.N. and Reinsch, C., "Balancing a Matrix for Calculation of Eigenvalues and 
Eigenvectors”, Numerical Math, 13, pp. 293-304 (1969). 

[27] Bavely, C.A. and Stewart, G.W., "An Algorithm for Computing Reducing Subspaces 
by Block Diagonalization” . SIAM J. Numerical Analysis. Vol 16, No. 2, April 1979. pp 
359-367. 

[28] Vivian, H.C. , Blaire, P.E., Eldred, D.B., Fleischer, G.E., Ih, C.-H.C., Nerheim, N.M., 
Scheid, R.E. and Wen, J.T., Flexible Structure Control Laboratory Development and 
Technology Demonstration , JPL Publication 88-29, October 1, 1987. 


106 


[29] Westreich, D. , ” A Practical Method for Computing the Exponential of a Matrix and its 
Integral”, Communications in Applied Numerical Methods , Vol 6, pp. 375-380 (1990). 

[30] Wie, B. and Bernstein, D.,”A Benchmark Problem for Robust Control Design,” Pro- 
ceedings of the 1990 American Control Conference , May 23-25, 1990. 

[31] Ge, Y., Colllins, E.G., Watson, L.T., and Doris, L.D., ’’Minimal Parameter Homotopies 
for the L 2 Optimal Model Order Reduction Problem”, submitted for publication in 
IEEE Transactions on Automatic Control. 

[32] Kleinman, D.L., Fortmann, T., and Athans, M., “On the Design of Linear Systems 
with Piecewise-Constant Feedback Gains,” IEEE Transactions on Automatic Control, 
Vol. AC-13, pp. 354-361, 1968. 

[33] Chen, B.M. , Saberi, A., Sannuti, P. and Shamash, Y., “Construction and Parameter- 
ization of All Static and Dynamic // 2 -Opt,iinal State Feedback Solutions,” Proceedings 
of the 31st IEEE Conference on Decision and Control, December 1992. Also, to appear 
in IEEE Transactions on Automatic Controls, March 1993. 

[34] Heffley, R.K., “A Compilation and Analysis of Helicopter Handling Qualities Data”, 
NASA Contractor Report 3145, Part II, August 1979. 

[35] “Aeronautical Design Standard” (ADS-33C), Handling Qualities Requirements for Mil- 
itary Rotorcraft , U.S. Army Aviation Systems Command Directorate for Engineering, 
St. Louis, MO, August 1989. 

[36] ElGhaoui, L., Carrier, A. and Bryson, A.E., Linear Quadratic Minimax Controllers, 
Journal of Guidance, Control, and Dynamics, Vol. 15, No. 4 (July-August 1992), pp 
954-961. 

[37] Khargonekar, P.P. and Rotea, M. A., “Mixed H 2 /H 00 Control: A Convex Optimization 
Approach,” IEEE Transactions on Automatic Control, Vol. 36, No. 7, pp. 824-837. 

[38] Chen, B.M. , Saberi, A. and Ly, U. “Simultaneous // 2///00 Optimal Control: The State- 
Feedback Case,” Proceedings of the 1992 AIAA Guidance, Navigation and Control 
Conference, August 10-12, 1992. 

[39] Doyle, J.C., and Stein, G., “Robustness with Observers”, IEEE Transactions on Au- 
tomatic Control, AC-24, pp. 607-611, 1979. 

[40] Ly, U., “H 2 - and //°° -Design Tools for Linear Time Invariant Systems”, Proceedings 
of the 3rd Annual Conference on Aerospace Computational Control, Volume I, Dec 
15,1989, JPL Pub. 89-45, Vol I. 

[41] Ly, U., VanSteenwyk, B., and Schomig, E., “Robust Control Design Using Parameter 
Optimization Techniques”, Control and Dynamic Systems, C.T. Leondes, ed., Aca- 
demic Press, 1993. 


107 


108 


Appendix A 


Preliminaries to the Hybrid 
Algorithm 

A.l Introduction 

It has been shown that the original gradient computation based on the diagonal form is fast, 
though inaccurate for defective systems. While the robust form described in Section 3.3 is 
accurate under most circumstances, it is also slower by an estimated factor of 8 or 9. The 
next question is whether one can apply the faster algorithm to the non-defective part of 
the matrix and the slower one to the defective part, thereby resulting in a fast and robust 
overall numerical scheme. Until such a scheme is actually implemented and tested, one can 
only speculate. There is a clear direction, however, as to how to structure that procedure. 

In order to have control over the separation of the diagonalizable part and the defective 
part one must abandon the simple eigenvalue-eigenvector decomposition for a combination 
of the Schur decomposition and a procedure to separate the set of defective eigenvalues from 
the non-defective ones. This custom routine can then feed the actual hybrid computation 
algorithm with the appropriate parts of the input matrix. 

Once a matrix is in Schur form, there will be a block diagonal (consisting of individual 
real eigenvectors and 2x2 blocks for complex pairs) along with an upper triangular part. 
The 2x2 blocks are not yet in a — u form. Moving non-defective eigenvalues apart from the 
defective ones usually requires zeroing the upper-triangular interaction between these roots. 
In fact, the very definition of a defective root depends on the ability to zero this cross-term. 


A. 2 Converting a 2x2 block to o — u form. 


Conversion of a general 2x2 diagonal block to the form: 
and a scale conversion, 


depends on a rotation 
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where s = sin(<p ) and c = cos(4>). The result of the rotation is to bring the two diagonal 
elements to equality. The rotated matrix becomes 

c 2 a 4- sc(f + 6) + s 2 g -s 2 f 4- sc(g - a) + c 2 b 
c 2 f 4 sc(g — a) — s 2 b s 2 a — sc(f 4 b) 4 c 2 g 


Thus the condition on the diagonal elements results in 


c 2 a + sc(f + b) + s 2 g = 
2 sc(f + 6) = 

sin(2(f ) ) = 2sc = 
cos(2(f)) = <? — s 2 — 


s 2 a — sc(f 4 b) 4 c 2 g or: 
(c 2 - s 2 )(g - a) 

(ff-a) 


A 


if+3 

A 


where 


A = v'a + ^ + Cs-a) 2 


Given values for sin(tfi) and cos(4>), it remains to find the scaling terms that will bring the 
off diagonal terms to being the opposite of one another (since the scaling will not affect the 
diagonal terms). The result is: 


or 


d2 ( (f-b) _ AN 

di\ 2 2 ) 


di / (/ -6) AN 
d 2 \ 2 2) 


da 

d i 


\ 


il , A 

2 '2 

A' 

2 2 


Of course the argument under the square root must be positive, which can be written as 
(/ — b) 2 > A 2 . 


A. 3 Reduction of a 2x2 Block to Upper Triangular Form 

The transformation from an upper Hessenburg form to the real Schur form for a matrix 
depends on the QR iteration. This iteration has some difficulty for a matrix with defective 
degenerate eigenvalues. Thus, there can be nonzero sub-diagonal elements not associated 
with complex eigenvalues. It is possible for one to see nonzero subdiagonal elements for 
a whole Jordan block, but usually one sees only a scattering of 2x2 diagonal blocks. It is 
sometimes necessary to check each 2x2 block to see if it is non-converged or if it really is a 
complex pair. Zeroing a sub-diagonal element can be done with a Givens rotation 


c s 


a 

b ' 


c —s 

—s c 


. / 

9 


s c 


resulting in 


c 2 a + sc(f + b) + s 2 g —s 2 f + sc(g — a) 4- c 2 b 


' a' (/' 

c 2 f + sc(g - a) - s 2 b s 2 a - sc( f 4 b) + c 2 g 


° g' 
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Thus, c 2 / + sc(g — a) — s 2 b = 0. This can be manipulated into 

cos{2<p) + sin(2<j>) = (~2^) ' 

To solve the above, we can consider it like a linear differential equation- there is a homo- 
geneous part and a particular part. We can write: 

cos{ 20) = Ki{a-g) + K 2 (f + b) 
sin( 2<j>) — K\(f + b) + K 2 (g - a) 

with 

K = tzl 

2 (/ + b) 2 + (g - a) 2 

[(g - a) 2 + 46/] 1/2 
{(f + b) 2 + (g-an 

K 2 is associated with the “particular” portion and scales it to be equal to the term 
K\ normalizes cos 2 {2<p) + sin 2 ( 2<j>) = 1 . Note that one has the condition: (g — a) 2 + 4b f > 0 
which complements the condition facing the transformation converting 2x2 blocks into the 
a — u form. In other words, if this condition is not met, it is a complex eigenvalue pair. 


A. 4 Reducing the Schur Matrix to the Diagonal Form 


One can subdivide the problem of zeroing the upper triangular portion of a Schur matrix 
into 4 cases that arise when one tries to zero the cross-term between roots. All of these cases 
are just specifications of the general form specified in [27]. When a matrix has been Schur 
decomposed there will, in the most general form, be an upper triangular partition along 
with a diagonal of real roots and 2x2 complex pairs. A normal matrix will be diagonalized 
by a Schur decomposition— the Frobenius norm of the strictly upper triangular partition 
of the matrix is known as a matrix’s departure from normality (see [23]). To diagonalize the 
upper triangular matrix we have to use nonunitary similarity transformations. The most 
convenient form is 


R = 


I P 

0 I ’ 


R~ l 


I -P 
0 I ’ 


Suppose that one were to zero the upper triangular part between two diagonal submatrices 
with this transformation: 


I -p ' 
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In general, one has the Lyapunov equation: T — PA 2 T AiP = 0 to solve. This equation has 
a solution provided that there are no eigenvalues in A t equal to those in A 2 . However, in 
general, the presence of degenerate eigenvalues, whether defective or not, will make deter- 
mination of this transformation impossible. Not only that, but in numerical computation 
defectiveness is more of a condition to be approached in a continuous manner rather than an 
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on-off setting. A measure of this more continuous form of defectiveness is in the condition 
of the transformation (diagonalizing) matrix. The condition of the overall transformation 
matrix will be helped by limiting the condition allowed for these individual non-unitary 
transformations. By declaring a minimum allowed condition, one thus detects, in a sense, 
which eigenvalues are part of a diagonalizable submatrix and which are not. (Because the 
inverse is explicitly known, this condition is easy to estimate.) 

It may actually be that a pair of distinct (but close) eigenvalues will get declared a 
degenerate pair, but this will only impose a slight time penalty on the eventual calculations. 
It can be said that the resolution on the eigenvalues was not sufficient to consider them as 
distinct. 

Fortunately, non-defective degenerate eigenvalues form a sort of normal submatrix so 
that their cross-terms are already zero as a result of the Schur decomposition — zero being 
the effective zero for a given machine precision, the oo — norm of the matrix, and the 
machine zero given the algorithm for the subdiagonal elements. It is possible to detect this 
relative zero condition and avoid problems in solving for a transformation matrix — by not 
having to solve for a transformation matrix at all. 

The end result will be the diagonalization of the non-defective eigenvalues, the con- 
tinuance of an upper triangular block for the defective eigenvalues, and a transformation 
that will remain well-conditioned. What becomes necessary, however, to keep the problem 
tractable is to break the cross-term zeroing into 4 cases: 2 real roots; a real root and a 
complex pair; a complex pair and a real root; and finally, between two complex pairs. 

A. 5 Reducing the Upper Triangular Part Between 2 Real 
Roots 

Reducing the upper triangular part between 2 real roots is the first (and easiest) case. One 
literally has 


1 -p 
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0 1 

L 0 
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resulting in p — . One does not need to solve for p if t \2 is relative zero for the matrix 

(thus, it would not matter if Ai = A 2 ). 

A. 6 Upper Triangular Part Between a Real Root and a 
Complex Pair 

This case is elucidated separately mostly as a means of showing notation. The submatrices 
of interest would be 
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Zeroing t\\ and £12 results in the following system of equations 


An — An 

A21 


P11 


tn 

A12 

A22 - An 


P12 


t\2 


which should always be solveable. In any case, one could always check for a zero solution 
here as well. 
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A. 7 Upper Triangular Part Between a Complex Pair and a 

Real Root 

This case is distinct from the above due to a slight difference in equations and notation 
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Zeroing in and i 2 i results in the following system of equations 
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A. 8 Upper Triangular Block Between 2 Complex Pairs 

Here one solves for the 2x2 block that zeros the upper triangular portion between a set of 
complex pairs (not necessarily in a — uj form) 
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The accompanying system of equations is 
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A. 9 Using These Techniques for Zeroing Blocks 

The emphasis here is to make sure that one can separate non-defective eigenvalues (and as- 
sociated eigenvectors) from a blocks containing defective ones. While physically separating 
individual eigenvalues (or complex pairs) from each other and from defective degenerate 
blocks while zeroing the upper triangular part between the two may sound appealing, it 
is easier from the bookkeeping standpoint to keep the diagonal terms fixed , then organize 
them at the end. 

Thus one can proceed through the upper triangular part of the matrix in a sequential 
order, using some sort of indexing array to encode the solitary eigenvalues, the complex 
pairs, and those eigenvalues belonging to a block of defective degenerate eigenvalues. At 
each upper triangular element (or block if this element associates one or more complex 
pairs), one can attempt to zero the members using the methods previously indicated. If a 
defective degenerate condition is indicated, the association with the eigenvalues will indicate 
the appropriate block to assign the eigenvalues. 

Those eigenvalues whose associated upper triangular parts have been zeroed can be 
moved past one another in a matrix by a simple shift transformation. The rest will collect 
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into blocks when they encounter eigenvalues they cannot shift past. These subblocks will 
be upper triangular and there will be no need to put them into Jordan form. 


Schur Form: Decomposed: 
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Appendix B 

Cost and Gradient Computation 
Using The Hybrid Algorithm 


B.l Introduction 

The current two methods within SANDY represent two extremes in the handling of matrix 
functions: the faster original “diagonal” method decomposed the argument matrices into 
eigenvalues and eigenvectors, then applied scalar functions to the diagonal terms; while the 
nominallly more robust method used an exponential series to compute a function of the 
whole matrix. 

While the diagonal method depends on successfully determining an eigenvalue-eigenvector 
decomposition, the robust form also has a weakness — precisely where the diagonal method 
is best suited. If one has scales in a matrix that are wildly varying, the robust form is 
more prone to error, even to the point where a defective-degeneracy does not look so bad. 
Consider a root matrix 
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matrix is transformed by 

the following matrix: 





" 0.2113 

0.4524 

0.6538 

0.7469 

0.1167 

0.2260 ‘ 



0.0824 

0.8075 

0.4899 

0.0378 

0.6250 

0.8159 

'T 

i 

0.7599 

0.4832 

0.7741 

0.4237 

0.5510 

0.2284 

1 


0.0087 

0.6135 

0.9626 

0.2613 

0.3550 

0.8553 



0.8096 

0.2749 

0.9933 

0.2403 

0.4943 

0.0621 



0.8474 

0.8807 

0.8360 

0.3405 

0.0365 

0.7075 


The condition of this matrix is 51. This should not introduce any complications in this 
comparison of methods. Because fundamentally all calculations depend on taking an expo- 
nentiation a good test case can be founded on the comparison of exponentiation calculations. 
Matrix exponentiation of the component submatrices of A and transforming by TAT -1 gives 
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the most accurate result, 


1.0E + 8 x 


' 0.4628 
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. 0.5702 
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0.8333 
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-0.3651 . 


Suppose one transforms A and then exponentiates, the difference with the previous calcu- 
lation would be 


’ -158.9406 
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56.7881 

-83.6180 

141.6097 

-81.9972 

-155.5265 

24.5134 

58.7931 


This error is still only on the order of one part in 10 6 . An exponentiation arrived at by the 
diagonal method does not have as much error though the matrix is defective degenerate, 


' -19.0856 

23.9609 

-13.3939 

-24.0870 

4.2366 

10.6407 ' 

-6.1430 

5.6210 

0.5094 

-6.3906 

-2.4559 

2.9409 

-15.4594 

17.7026 

-6.9271 

-18.3993 

0.3217 

8.2278 

-15.4252 

16.3490 

-3.9463 

-17.4740 

-2.0148 

7.9009 

-15.4366 

16.2809 

-3.7615 

-17.4420 

-2.1720 

7.8944 

. -15.2599 

17.0281 

-5.8801 

-17.8393 

-0.4249 

8.0123 


The condition of the eigenvector matrix is 9.9927 E + 4 due to the defectiveness, thus the 
calculation is probably not too impaired. 

Thus, it is hoped that, in addition to an increase in calculation speed, this algorithm 
will be made more complete. 


B.2 The X Integral 


The simpler of these integrals is X{t) = e Ar Be Cr dr. The matrices A and C can be de- 
composed along the lines of an eigenvalue — eigenvector decomposition, but only for the non- 
degenerate eigenvalues. The degenerate eigenvalues are decoupled from the non-degenerate 
eigenvalues, but are otherwise left in a non-diagonalized matrix. The non-diagonalized part 
of the matrix will consist of decoupled blocks of degenerate sets of eigenvalues. This decom- 
position is based on a Schur decomposition with a selective eigenvalue shift and a decoupling 
algorithm. We have 


A = V A 


' A a 

0 

V7 1 and C = V C 

Ac 

0 

0 

Wa 


0 

Wc 


Vc 1 


with Aa and Ac both being diagonal. For the exponential 

V ^‘4 = VAexp 


e At = exp < Va 


A A 0 
0 W A 


{ 


A a 0 
0 W A 


t > v\ 


-1 
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This results in the following form for X 


X{t ) = i e Ar Be Cr 

Jo 


dr 



A a 0 
0 W A 



Ac 0 
0 W c 


rjdr j^ 1 (B.l) 


where B = V^ l BVc- The exponentials can be expanded and X can be separated into parts 
involving only the non-defective eigenvalues, the defective degenerate part, and combina- 
tions of both. For 


Bn B n 
B 21 B 22 


(B.2) 


we have 


X(t) = V A 


j} e A * T B n e AcT dr J* e A * T B l2 e WcT dr 
&e WAT B 21 e AcT dr / 0 ‘ e w * T B-ne WcT dr 


(B.3) 


As a result one can use the algorithm used in [18] for the upper left partition (usually all or 
most of the matrix), and thus preserve speed. The remaining portions involve various uses 
of the exponential to represent the integral. 

Consider the similar forms for the integrals of the upper right and lower left blocks. 
These involve a combination of a diagonal portion and an irreducible portion (typically of 
dimension much smaller than the original matrix). For example, one can rewrite the upper 
right integral in the following way 


with 


*w(4) = r 

Jo 


e AAT B l2 e WcT dr 


(B.4) 


X kw (t) = the k th row of *12 


= f 


e XAkT b kW e WcT dr 


(B.5) 

(B.6) 


For when \ A k is a scalar, b k w is the k th row of $ 12 - The number of columns in b k w 
corresponds to the number of columns of the degenerate matrix Wc- Because these matrix 
integrals can be broken down row by row, in the case of the upper right block, or column by 
column, in the case of the lower left block, there is a considerable time savings. Computation 
of the matrix exponential corresponding to the integral Jq e Ar Be CT dr for a general A , B, 
and C is on the order of ( dim(A ) + dim(C )) 3 floating point operations. Breaking down this 
integral by rows (or columns) results in the advantage of evaluating several small integrals. 
This special integral form allows us to write: 


X kw ( t)= [\ kW e( Wc+XAkI > dr 

Jo 

An algorithm for evaluating the matrix exponential integral will be presented in section B.4. 

For A Ak being a 2x2 complex root pair block, the situation is a bit more complicated. 
Since 


a 

u 

= 6* 

COS U) 

sincj 

—w 

a 


— sin a; 

COS U) 
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one can expand the integral (B.6): 


4;(i)Ml] = / e° AkT cos(u)AkT) b k we WcT dr + [ e° AkT k r) b k +\we WcT dr 

Jo Jo 

Xkw{t)[col2] = f e aAkT cos(uAkT)bk+\we WcT dr — [ e tr/,fcT sin^^fcr) bk\ve WcT dr 
Jo Jo 

Again, with e aAkT , cos u^r, and sinter as scalars we can re-arrange the above integrals 
into forms depending on two basic terms: 


and 


f 

f 

Jo 


e (W C +tTAkI)T COS (^ U J AkT 'j d-j- 


e (Wc+*Aki)T s [ n ( UAkT ) dr 


An algorithm for computing these terms will be presented in section B.5. 

The lower right portion (containing only the degenerate eigenvalue matrices) will have 
to be evaluated with the pure exponential form of the integral evaluation. The size of the 
degenerate matrix blocks should not be very large relative to the number of non-degenerate 
eigenvalues. Not only will this keep the overall computation from slowing down, but the 
eigenvalues shall nominally be of the same relative magnitude. 


B.3 The M Integral 


This integral yields a somewhat more difficult set of cases of integrals to evaluate, as we 
shall see. 

The integral M(t) = Jo Jo e A< ^ v ~ s ' > Be Cv De Es ds dv becomes, when A, C, and E are 
partitioned into non-defective and defective parts. 


M{t) 


rJ’f 


6 Aa(v-s) 

0 

B 

' e A cv 

0 

0 

e W A {v- s ) 

0 

e w <*> 


■ e A £ s Q 

0 e W£S 


ds dv T e 1 


f l [ V e AA ( v ~ s ) B n e Acv 
Jo Jo 

| ^ [ V e AA(v ~ 3) B n e Acv 

Jo Jo 

T>ne A£S ds dv 

| V 12 e WES dsdv 

+ [ [ V e A ^ v ~ s) B 12 e Wcv 
Jo Jo 

| + T f V e AA ^ v - s) B 12 e Wcv ; 
Jo Jo 

T> 2 \e A£S ds dv 

1 V 22 e WES ds dv 

t f V e WA{v ~ s) B 2 ie AcV 
Jo Jo 

| 1* I” e WA <- v ~ s ) B 2 ie AcV 

Jo Jo 

T>ne AES dsdv 

| V 12 e WES dsdv 

i 

Jo Jo 

+ f f V e w ^-») B 22 e WcV 
Jo Jo 

V 21 e A£S dsdv 

| V 22 e WES dsdv 
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where 


TX l BT c 

T^DTe 


= B = 
= V = 


Bn &\2 
B21 B22 
T>n V \2 
T> 2 \ T >22 


The first part of the M\\ summand can be evaluated in the standard fast way, so it 
would be desirable to be able to evaluate the second integral as quickly. One can solve 
analytically for the integral over s. Let’s assume initially that the eigenvalue is real. We 
have 


f V e Xi{v - s) Bn2e Wcv V2ije x i s dsdv 
Jo Jo 


Jo <M ' : 0 T e<Ai_Al) ' ds ) Bine Wc ' 1 V 2 \, dv 

i 


rt p^jv _ pXiV 

-B il 2 e WcV V 2 u dv 


A j A 1 


g * 12 f* f e (W c +\jI)v 
A j J 0 ' 


_ e {W c +\il)v^ dy V2i . 


The end result is just the matrix exponential integral again. For a complex pair, one has 
to deal with a 2 x 2 block (looking just at the integrand) 


e <Ti(f S) CQSUIi ( V _ S ) g il2 + e CT<(u-3) — S ) Bi+1,12 

e < 7 i(ti-s) cosu) .( v _ s ) Bi+i t i2 — e CTi ( v ~ s ' ) smu>i(v - s) Bn2 


e W c v 


cosujjS | T>2ij+ie <TiS cos u)jS 
— T>2ij+\e <T i s sin UjS | V 2 ije CTiS cosu-jS 

This becomes 


e aiV e (crj-<Ti)s 



cos Ui(v - s ) B l i2e WcV T>2ij 

| [coso;i(u - s ) Bii 2 e WcV T>2i t j+i 

+ sinui(v - 5 ) B l +i f i 2 e WcV T>2ij] 

| + sin Ui(v - s) Bi + i t i2e WcV T>2ij + i 


•COS UjS 

| -COS LJjS 

— 

cosuJi(v - s) B i i2e Wcv T>2i l j+i 

| + [cosWt(u - s) Bn2e WcV T>2\j 

+ sinu;i(u - s) B i+ i i i2e WcV T>2i tj +i 

| +sm(jji(v - s)Bi + i'i 2 e WcV T> 2 ij 


-sinujjS 

| • sinews 


)cosuJi(v - s)B i +i i i 2 e Wcv T>2ij \ 

- sin LJi(v - s ) Bii2e WcV T>2ij} | 

•COSWjS I 

- j^cos u)i(v - s) Bi+i'i2e Wcv T>2\j+i | 

-smui(v - s)Bn 2 e WcV V2\j + i \ 

■sinWjS | 


[coswi(u - s ) B i+ i,i2e Wcv T>2ij+i 
-sina)i(v - s) Bn2e WcV T>2ij+\ 

■ COS U)jS 

[coswi(u - s ) B i+ i t i 2 e Wcv T>2ij 
— sino)i(v - s) Bine WcV T>2ij 
• sin UjS 
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The integrals over s will boil down to these terms 

Jq cos (uj - Ui)s ds £' cos (uj + u^s ds 

/o e^J ~ a ^ s sin(wj — cj 1 )s ds Jq e^ a i~ <T ' >s sin(wj + w,)s ds 
Substitution of these back in for our 2x2 matrix will yield a series of integrals of the form 


and 


/ 

f 

Jo 


e (W c +<r'I)v CQg u) r y ^ 


e (W c +cr'I)v s j na / w d v 


The amount of algebra to get this matrix into a series of integrals on v is substantial. 

Being able to re-express the integral over s in analytic form is not only sufficient for 
simplifying the double integral but necessary. Changing the order of integration will in- 
troduce numerically unstable terms that can be handled only by doing the whole double 
integral numerically in the current “robust” M calculation. Thus one needs to have only 
simple eigenvalues involved in the s integral. This happens only on the integrals in the M.\\ 
partition. The rest of the integrals will have to use the current robust form in some manner. 

One could be given the impression that not much work and time is being saved. However, 
a close examination of what happens when one does these other partitions with the robust 
algorithm one will get a flavor for the amount of time saved. 

Consider the sum for M 22 . With the following 

/* f e WA{ ~ v - s) B 2 xe KcV V n e WES dsdv = V /* /" e WA(v - 3) B 2U e XciiV T> il2 e WES ds dv 
Jo Jo “Jo Jo 

Considering that the overall cost of a direct matrix exponentiation is of 0(n 3 ) or more, 

there is a net win. The fact that A c u may represent a complex pair in the form a W 

—U (J 

is of no consequence as the summands will be evaluated by a direct matrix exponentiation. 


B.4 Pade Series for Matrix Exponential Integral 


In the hybrid formulation, one of the terms boils down to finding the integral of the expo- 
nential of a matrix / 0 £ e Ar dr. It was suggested [23] to compute this integral by taking the 
following exponential: 


exp 


0 I 

0 A 


I foe Ar dr 
0 e At 


This technique was disputed in [29] , where it was determined that a direct Taylor Series for 
this integral was considered a better method. The basis for this was that the larger overall 
matrix contributed not only to a lengthly computation time, but also a loss of accuracy. A 
third alternative is presented here in the use of a Pade Series for computing the integral. 

One could call this series the generating function for the integral, but for the fact that it 
is simpler to assume that one is integrating over a unit interval. Since the matrix contains 
all the scaling information in a typical function call, it is a safe assumption. We know that: 



e 47 dr = 


A~ l e At 1 = A~ l \e A — il 
o 1 J 
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and 


e 


At 



E,=o D t (Aty 


SO 


f 


e Ar dr 


-x YSUFi-DQW 

E,“ 0 Di(Aty 


The most commonly used case of the Pade Series for the exponential is a diagonal sequence 
(order of the numerator is equal to that of the denominator). For this, a nice simplification 
in the above equation occurs since, for i even, Ni = Di ) and for % odd, TV* = — ZV Thus, we 
can write the Pade Series for the integral generating function as: 



dr 


tP ,N M {z) = 2t 


g^M-i+2 

E,= 0 o^Aty 


with 


Ni 

Di 


(2N -i)\N\ 

(2 N)\i\(N-i)\ 
t (2 N-i)lNl 
{ ; (2 N)\i\(N-i)\ 


With t being 1, we have the integral being equal to this particular Pade Series. 

Nominally one has to scale the input matrix so that its oo-norm is less than 1/2. The 
rescaling of the result is somewhat different than that of the matrix exponential. The matrix 
itself doubles in scale, yet the integration interval does not change 

r\ 1 r 2 

[ e 2At dt = - f e Ar dr 

Jo 2 Jo 

- Wo eATdr+ J! eMdT . 

The term e A is easy to generate from the integral at each step, i.e., multiply fg e At dt by 
the matrix A and add I. 


B.5 Pade Series for Integral of Matrix Exponential and Si- 
nusoid 

Another special form to handle from the hybrid formulation is 

ft ft 

/ e Ar sin uit dr or / e Ar cos ojt dr. 

Jo Jo 


Using the relationship 



3 

b 

IT 

COSCJ 

sincj 

exp 


— e 



—U <7 


— sin u 

COSCJ 
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Thus we can write 

W + al ul _ e (w+«r/) o I cosu) Ismu 

GX ^ — ul W + a I 0 e (W+ai) —/sin u /cos u 

__ e( w+cr ^ cos u e^^^sinw 
—e( w+aI ) sinw e ( w + aJ ) C os u 

It would seem that we are doing twice as much work as we need be. Since the augmented 
2 n x 2n matrix has several symmetries, one can do specialized sets of multiplications to 
exploit these symmetries and avoid use of more matrix storage than necessary. In the series 
calculation, one will need to take squares, etc of the argument matrix 

'W + al ul "| 2 T (W + oI) 2 -u 2 I ucrl 

-ul W + (t1\ ~[ -uol {W + al) 2 — u 2 I 

One will have only 2 distinct submatrices within this matrix. Thus one need not do explicit 
multiplies. 

B.6 Scaling 

Because all matrix blocks will have their eigenvalues on the diagonal and no sub-diagonal 
elements, a single scaling parameter characterizes the magnitude of the resulting exponential 
of the block. For example, in the block 



where Ai through A n are the same value. Because 

e wt _ e (A i+w')t _ e *t e w't 

one can start with an integral like 

[ l e WAT Be WcT dr 
Jo 

and “balance” the two matrices into having the same diagonal value. Note that the smaller 
modulus diagonal terms will probably reduce the number of scaling steps and/or series 
length in the matrix exponential calculation, thus saving some time. 
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For the calculation for M, the balancing of the matrix blocks is a more complicated job, 
due to the existance of two separate integrals involving separate matrices 



e A{v ~ a) Be Cv De Es dsdv. 


However, because the procedure is simple, the attempt is worth it no matter how little 
results one gets. 
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Appendix C 

An Algorithm for SCB 


The following sections show different steps in the computation of the SCB transformation 
discussed in Chapter 5. It is hoped that by proceeding from the simple cases and ending 
with the more general ones, a better understanding of the algorithm can be achieved. 


C.l S.C.B. for Non-Strictly Proper SISO Systems 


In a SISO system, the presence of a direct feedthrough term essentially eliminates most of 
the work needed in computing the special coordinate basis. In this case, for the system 


{ x — Ax + Bu 
y — Cx 4- Du 


(C.l) 


where u and y are scalar. With D ^ 0, all the states are basically “internal” states -the 
input can go to the output directly without necessarily passing through an integrator. One 
can rewrite the state equations by eliminating u in terms of x and y as follows, 

x = {A- BD~ 1 C)x + BD~ l y 


Clearly there is no longer a direct input term from u to any of the system states once the 
substitution has been made. Note that in this case, the system invariant zeros are simply 
the eigenvalues of the matrix (A — BD~ 1 C). 


C.2 S.C.B. for Strictly Proper SISO Systems 


Conceptually, this is the next easiest case. With no direct feedthrough term, the task is to 
determine an output state with a direct feedthrough term through differentiation. Because 
there is only one input and one output, only “output” states can have a direct input term. 
Furthermore, the output states or its derivatives will always have a direct input term, 
indicating that a SISO system is always left and right-invertible. We begin with the system 


{ x = Ax + Bu 

y — Cx 


(C.2) 


where u and y are scalar. Suppose that one transforms the states using an orthogonal 
transformation V c so that one of the states is the output, 


Cq = [c 0 ... 0 ] = CVc\ Aq = VqAVc — 


Aq yy Aoyx 

Aoxy A^xx 
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and 


Bo = VgB = 


bOy Q 

bO Xl 


resulting in the state equations 


b 0 X 


(n-l) 


Boy 

Bq x 


'y o' 


A()yy 

A()yx 

'y o' 

i 

B 0 y 

±0. 


. Aoxy 

Aq xx _ 

x 0 _ 


. Bq x 


In essence, we have taken one derivative of the output and this also implies that the system 
has at least one infinite zero. 

If the term Bq v (a scalar in a SISO system) is zero, then we continue to generate another 
output state by differentiating yo. Namely, 


j/0 AoyyijO + A()yxXO 

Making the appropriate substitution for Xq, we obtain 


Vo = Aoyyi/o 4- Aoy X [AqxxXq + A 0xy y 0 + Bq x u] 


Because yo is now a state variable, we can write the state equations for y 0 , yo as 


V i 


0 1 
AoyX Aoxy Ao yy 


Vi + 


0 

AoyxAo X x 


Xo + 


0 

AoyxBox 


u 


where 



y o 
V 0 


Note that in SISO systems, yi always contains the term with the highest derivative, and 
2/i covers the rest of the output states. In addition, one of the internal states xq can be 
rewritten in terms of yo, yo and the remaining states in Xq using the equation 


J/0 ^OyyyO "h ^Oyx^O (C.3) 

Identifying which state among x Q to be replaced by y 0 depends upon the structure of 
Aoyx — it is prudent to choose the state corresponding to element in Aoy X having the largest 
magnitude. Suppose that the j th element of the state x 0 is chosen. Let’s denote this state 
by xq and the rest of the internal states in xo by xo with 


%o\ 

202 


xo = 


X 0(j+1 ) 


L ^0(n— 1) J 
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Partitioning the state equation for xq into separate equations for x 0 and Xq, we have 


Xq — AoxxXo “I" Ao x £Xq p AQxyVO P B() x a (C.4) 

Xo ~ AqxxXo p AqxxXo p AoxyVO P B()x a (0.5) 

Then using equation (C.3), we can express xo in terms of the output states y\ and the 
remaining states x\ as follows, 


2/0 A()yy JJO A()yxXo — Aq 7 j X Xq , 


or 


where we define 


20 = \ l x ([- A Oyy I] Vl ~ AoyxXx) 


f 


y i 


y o 

V1 

— 


— 


in 


yo 

. Xi 

= Xo 




(C.6) 


So far we have replaced two of the system states by output states. The next part of the 
algorithm can be done recursively. 

First, we can re-write the system model in terms of the output states y\ and the re- 
maining states x\. The resulting state equations are 


where 


Alyy 

A\yx 

A l xy 

A hi 

B\y 


B\ x 


\ y\ ~ A\ yy y\ + A\yxX\ + B\yU 

| Xi = A\ xy y\ + A\ xx xi + Bi x u 


0 1 
Aoyx Aq x y Aftyy 

0 

A()y x A()xx 


P 


0 

AlOyxA 0: 


0 

A[}yx A()xx 


Aojx [~Aoyy I] 
Aoy X Aoyx 


[^Oxy P Aqxx ^Oyi [ A(iyy /]j 

Aqxx Aq xx A(yy X /!()?; i: 

0 

-4()yx Box 

Bq± 


(C.7) 


is essential for starting the next iteration. There is no point 

of taking the derivative of yi in the formulation of new output states, as they have no direct 
interaction with the X\ states. Recall from the first step that we have y\ = yo- Subsequent 
iterations would involve repeated differentiation of y\. Finally, at the /c £/l -iteration the 
output state vector yk is given by 


Partitioning of y\ into 


y i 
in 


where the output yk represents the k th derivative of the output yo. 
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At the k th iteration, we have 


f Vk 
\ x k 

where 


A kyyVk "h A kyx Xk -\- B k yU 
AkxyVk A AfcxxX/c "I" B kx U 


B k y 


A 


kyx 


A 


kyy 


0 

0 


0 

bky 


kx(n—k) 

[_ Akyx 

0 
0 

0 

. Akyy Akyy J 




(C.8) 


(C.9) 


(C.10) 


(C.ll) 


If the term B ky has some nonzero entry (i.e., b ky ^ 0), then the SCB transformation is 
complete. Otherwise, one must continue to take derivatives of the output state y k as long 
as there are nonzero entries in the term Akyx- Note that both the terms B ky and Ak yx 
cannot be zero simultaneously in a SISO system since it is always right-invertible. 

When the term B ky is nonzero, then the k th derivative of yo has a nonzero control input 
term. The control input u can be substituted out in terms of the output state, its derivatives 
and the remaining internal states x k as 


a — 1 /bky ilk AkyyitV Ak yy y k 


This can also be accomplished by performing the following state transformation 


where 

and 


Vf 


I 

0 ' 


Vk 

X a 


-P 

I 


Xk 


P = Bk x {Bl y B ky )- l Bl y = B kx [00 • • • 0 1 /b ky ] , 


Aff Afa 


I 0 

1 

<■¥ 

$ 

1 

7 0 

A a f A aa 


-P I 

1^ Afcxy A^kxX J 

P I 


After the transformation, the u input is removed from the state equations associated with 
the internal states x k and we arrive at the final SCB form, 


Vf 

X a 


Aff 
. A a / 



Vi 

Xa. 


+ 



U 


The system invariant zeros are simply the eigenvalues of the submatrix A aa . Clearly, some 
criterion must be established in testing whether the term bky at each iteration is zero or 
near zero. Often, a pre-specified level of tolerance must be given to the SCB algorithm for 
this singularity test. 
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C.3 S.C.B. for Multi-Input and Multi-Output Systems 


Let’s consider the following system 



Ax + Bu 
Cx + Du 


(C.12) 


where we assume without loss of generalities that the matrices [B T , D T ] T and [C, D] have 
full rank. 

Analogous to the non-strictly proper SISO case, we apply the singular value decompo- 
sition to transform the inputs and outputs so that the direct feedthrough term takes on the 
following form 


( 7 1 0 0 

0 02 o 0 


D = Ud'EdVd — Ud 


0 

0 




0 


vS 



0 


where Ud and Vo are orthogonal transformations and the singular values are given in a 
decreasing order { (J\ > d 2 > ■ ■ ■ > <i Tn > 0). The transformation Vp is applied to the input 
u while the transformation Up is to the output y. Further classification and partitioning of 
the x states will still be required, as well as for any remaining input and output terms that 
are not related to the direct feedthrough term. Note that the singular values can be zero 
for some k <m. With the above input-output transformations, we define 


and 


The output equation becomes 

'y 

.y 

Partitioning of { u , u } and { y , y } depends on the rank of D. The matrix £ /; represents the 
upper-left nonsingular matrix of And the matrix C and C are given by 

= uj>c 

We can express u in terms of the output y and the system states x from the equation 
y = Cx + t, D u as 

u = £p 1 [y-Cx\ (C.13) 


C 

c 


u 

u 


= v£u 


y 

y 


= ul~y 


C' 

c 


x + 


Sd 0 
0 0 


u 

u 
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The system given in equation (C.12) can be re-written as 


x = Ax + BVp 


u 

u 


or 


x = Ax + Bu + Bu 


Using equation (C.13), we have 


x = Ax + Bu + Ly 
< y = Cx+tiou 

y = Cx 

where 

A = A - B±~ D l C 
L = B ±~ d 1 


(C.14) 


First let’s consider the case where u comprises all of the inputs, hence the input u is non- 
existent. Furthermore, if the outputs y is also non-existent, then the SCB transformation 
is complete resulting in the following system, 


{ x — Ax + Ly 
y = Cx + £,du 


(C. 15) 


Clearly, the system is right- and left-invertible with no infinite zeros. Furthermore, the 
system invariant zeros are simply the eigenvalues of A. 

Let’s examine next the simple case where y comprises of all the outputs, then the 
remaining inputs u must be associated with the internal states x in the equation 


{ x ■= Ax + Bu 4- Ly 
y = Cx + Y, D u 


(C. 16) 


We can further distinguish the system states x into two subspaces identified with states x a 
and x c corresponding respectively to the uncontrollable and controllable subspaces of the 
inputs u, as discussed in Section 5.3. The resulting system is in the following SCB form 


x a 

X c 



Xq, 
X c . 


+ 


0 

B c 


u + 


L a 

L c 


y 


where the pair ( Acc , B c ) is controllable. Clearly, the system is right-invertible (but not 
left-invertible due to the presence of u) and has no infinite zeros. In this case, the system 
invariant zeros are simply the eigenvalues of A^. 

If on the other hand there are outputs y in equation (C.14), then they need to be 
converted into states in a manner similar to that shown in Section C.2 for SISO systems. 
The procedure will be described in more details in Section C.4. 

In the next section, we examine the general case where both u and y exist in the 
equation (C.14). 
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C.4 SCB Transformations for Decomposing the Remaining 
States 

Let’s consider the system given in equation (C.14) where the direct input u has already 
been eliminated through the output y. We now have the following system, 


x = Ax + Bu + Ly 
y = Cx 


(C.17) 


The task is to resolve the system into the 4 types of SCB states yj, yt,, x a , and x c . The 
first step follows conceptually the procedure in Section C.2 for SISO systems. We begin by 
performing the singular value decomposition on the matrix C, 

rT 


C = UcVVS 

= Uc[t c 0] Vq 
= [Co 0}V£ 


(C.18) 


where Uc and Vc are orthogonal matrices and S c is a nonsingular matrix of dimension 
p — mo which is actually equal to the rank of C. 

At the first iteration k = 0, we transform the states x such that some of the states are 
outputs with the following state transformation 


Vo 

Xq 


= VJx 


(C.19) 


Letting Uq = u, the state model in equation (C.17) becomes 

+ Vq Buq + Vq Ly 


y o 

X 0 


= VZAVc 


Vo 

X 0 


or 


l 

y 

= CV c 

o o 
5% H 






H* CC£. 
o o 

= 

A()yy A.Qyx 

_ A()xy A$xx . 

o o 

+ 

Boy 

. Bqx . 

"UQ + 

t-" 

0 o 
H « 

1 1 


(C.20) 


y 


= [Co 0] 


y o 

Xo 


(C.21) 


Next, we identify the output states that have direct interaction with uq by examining the 
term Bo y . Again using the singular value decomposition, we have 


Boy = U Boy ^B 0y VB 0y 

S B 0y 0 
0 0 




= u Bo v 

We then apply the following transformations to the inputs u and the outputs yo, 

'Vo' 


(C.22) 


and 


u BoyV o = 

VS»«o = 


Vo 

uo 

Uq 


(C.23) 


(C.24) 
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Clearly, 


Vo 

[Jo 


“ U B 0y A Oyy Ub 0v ^ + U s Qy A 0yx X 0 + U^ Qy B 0y V Boy 


VO 


Uq 

UO 


or 


and 


or 


V o 

£oJ 


Aoyy A 0 ^ 

A&yy A)yS J 


3/0 

SoJ 


+ 


Xo + 

^Boy O' 

Wo 

+ 

Loy 

- 

0 0 

.wo. 


- L 0y- 


®0 AoxyUgoy 


3/o 

IJo 


■^Oyz 


+ + Bq x V B( 


o» 


Wo 

WO 


Xo — [ Ao X y A 0x y j 


y o 

L£oJ 


+ ^ohXq 4- [ B 0x V Bov B 0x V Boy ] 


Wo 

Wo 


+ U T Boy L 0y y (C.25) 

(C.26) 
(C.27) 
(C.28) 


+ Lo x y 


+ Lo x y 


In the state equation for x 0 , we can eliminate the dependency on u 0 using the following 
equation 

“o ^B 0y [ 3/0 A()yy 3/0 *^0yy3/0 -^Oya:®0 f^Oyj/j 

This procedure can also be accomplished by invoking the following state transformation 


3/o 


' 7 0 0‘ 


y 0 

% 

= 

0 7 0 


3/0 

Xq 


1 

O 

ft, 

1 

- t 


Xo 


where 


P — BqxVbqJ^ B l 0y 

Let’s define the state matrices in the transformed system as 


-^Oyy -^0y5 Aoyx 
Jyy ‘^Oyy '^Oy:r 
A()xy ^Ozy Aqxx 



I 

0 

0 ' 


~ 

0 

7 

0 



-P 

0 

7 



Ao yy A(jy^ Aq£ x 
A ojjy A()yy Aog x 
Aoxy Aq x ^ Aqxx 



' 7 

0 

0 ' 


0 

7 

0 


P 

0 

7 


and 


1 

M> 

£ 

0 

L 0 y 


7 

0 

0 ' 


£ Soy 

0 

Loy 

0 

0 

L 0 y 

= 

0 

7 

0 
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0 

L 0y 

0 

Box 

Lqx 


-P 

0 

7 


. B 0x V b 0y 

Bq x Vb 0v 

Lqx 


[Co Co 0] = [do Co 0] 


7 0 0 
0 7 0 
P 0 7 


where [C 0 ,C 0 ] = C 0 U Boy . 

The resulting system is given by 


'y o' 

r 

h 

= 

-Xo. 

- 

y 

= [< 


1 0yy 


A&Qij -^Oyz 


Aoyy A)yy ^Oyz 
Aoxy Ao^y Aqxx 

~ Vo' 

t/o 

1x0. 



Jo' 

' 


~Vo 

+ 

. 

-Xo. 



E Bo„ 

0 

L 0 


0 

0 


Loy 

L Oy 


Box Lqx J 



"fio - 


UO 

. 

. V . 


(C.29) 
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Next, we identify in the state equation for y 0 those outputs that have direct influence from 
xo- The output vector y 0 is then partitioned accordingly into yo and yo- The partition is 
done based on the rank of A^. Similar to previous matrix decompositions, we perform 
the singular value decomposition on the matrix A^, 


A()yx ^OAyx^OAyx^oAyx 


Ez . 0 

n Oyx 

0 0 


(C.30) 


where U 0j and V 0 ^ s are orthogonal matrices. Let’s define the following state transfor- 
mation 

I“j=^0 (C.31) 

and the following output transformation 


V o 

L 2/o J 


(C.32) 


The output state equation for y is now partionned into the following output state equations 


j y o Aoyyfjo + A()yy 7/() + A()yyjif0 A ^OA^X^O A To yy 

Vo ~ Aoyyyo + Aoyyyo A Aoyyyo A Logy 


(C.33) 


where E 0 ,j« £ is a square and invertible submatrix of ^oAyx- 

By taking derivatives of the yo states, the algorithm converts system states xq into output 
states. This technique will be repeated on subsequent steps until the overlap indicated by 
the matrix Akyx is zero at some k th iteration. 

Note that if there is no direct feedthrough term Loy in the y 0 state equation, then one 
can simply take its derivative as follows, 

Vo = Aoyyy 0 4- AoyyXjQ + Aoyyf)o A Eq^^^O (C.34) 

thereby creating y 0 as a new se t output states replacing io. Of course, we need to 
substitute y 0 , yo > and £q from equations (C.29) and (C.33) in the above equation. The 
term y 0 is not substituted out because it is no longer part of the state equations — this 
equation is used to replace ’xo states with f/ 0 states. 

However, if the direct feedthrough term Loy is not zero, then the term Loyy would be 
present in the equation for y 0 in equation (C.34). It is impossible to create y since this 
term would include a time derivative of the control term u. However, instead of creating 
the states y 0 , one would consider a new output 

Vo = Vo ~ L oyV _ 

= Aoyyyo A Aoyyyo A Aoyyyo A ^oAyx^o 
Differentiating the above equation we have 

Vo ~ Aoyyyo + Aoyyyo A Aoyyy 0 + ^oAyx^o 
Note that we have a new output state equation from the above definition of y 0 

Vo = Vo + To yy (C.37) 


(C.35) 


(C.36) 
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and equation (C.36) becomes 

Vo — AoyyVo + Aoyy [y 0 + LoyV] + Aoyyy 0 + f^QAyx^O (C.38) 

Substituting the appropriate equations for y 0 , y 0 and x into equation (C.38), we obtain 
Vo = A 0yy [y 0 + L og y] 

+ A(iyy Aoyyi/Q + A^yy^Q 4" Aoyyj/O + AqQ^Q + Aq^^Q + T, B Qy UQ + L^yV 
+ Aoyy [Aoyyyo + AoyyVo + AoyyVo + ^OyV] 

4" ^-‘Aoijx [ A)xy Vo 4" A^^yo 4~ ^o^yVo 4“ Aq^^xq 4” A^^xo 4” Bq^Uo 4" ^o±y 
After simplification, we have 

Vo = ^oyyVo 4- A)yy yo 4- Aq^Vo + A Q ^y 0 4- A 0 ^k o + 4- B 0 ^uo + Bq^uo + L 0 ^y (C.39) 

Using equation (C.35), we can solve for xo in terms of y 0 , yo, yo, Vo and y. Namely, 

^0 = ^-‘oAyX [So 4" LoyV — Aoyyf/Q — AoyyljQ — Aoyyijo] (C.40) 

With equation (C.40), one can eliminate the dependency on from all the above state 
equations. The resulting system state model has the following form 

A)yy Aoyy 
0 0 


~Vo~ 


So 


Vo 

= 

Vo 


to. 



Aoyy A(jyy 


Aoyy A()vv 


Oyy 


0 

Aoyy 


1 

0 


yy . Oyy 
^Oiy A)iy 


Aoyy Aoyfi 


•^Oxy 


= [Co Co Co 0 0] 


•^Oxy 

V o 

yo 

yo 

Vo 

^o 


•^Oyx 
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0 

Aoy'x 
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'yo' 
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yo 
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E B o v 
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0 
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Loy 

L 0 y 

>_ 0 y ^Oy 
Box L q = j 


0 

0 

0 

Bn 


u 0 
i to 
L V 


(C.41) 

Let’s apply the following state transformation to eliminate the term B 0 § in the input matrix 
associated with the input Uo, 


where 
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(C.42) 


(C.43) 
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Applying the same transformation Tp to the system and output matrices, we obtain 


A 0yy 

A 0 yy 

Ao yy 

A - 

^0 yy 

Ao yx 


Aoyy 

Aoyy 

A 0 yy 

A 0yy 

Ao yx 
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— Ip 
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A 0 yy 

A 0 yy 
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0 

A 0~yy 

A - 

^ Oyy 

A 0 </y 

Ao yy 

A -- 

^0 yx 


A 0yy 

A 0yy 

A 0 yy 

A 0yy 

Aoyx 

- A 0£y 

A 0 ±y 

A Oxy 

A - 

^0 xy 

A 0±x - 


- A 0xy 

A 0±y 

Aq xy 

A 0xy 

AqM - 


[Co Co C 0 0 0] 

The new system is of the form 


[Co C 0 Co 0 0] Tp 

(C.44) 


'k' 
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- Ao^y A 0xy Aq xy 

><?> i 

y 

o 

o 

CJ 

II 

1 

o o o 

1 


A 0y± 


'yo' 


’ S B 0y 
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0 
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Uq 
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Vo 

+ 
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0 

L()y 


Uq 

A)$i 
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0 

A)y 

L 0i 


. V - 

A 0xx . 


Xq 


0 

A)1 

Loi- 



(C.45) 

The next iteration can be carried out in a similar fashion for the subsystem corresponding to 
the system states y 0 and Bq. For convenience, we re-define them as y i and X\ respectively, 


( 2/1 i° (C.46) 

[ X\ = xo 

and letting the inputs y\ and ui be 

‘ y 

yo 

y yo 
.yo 

ui = uq (C.48) 

At the iteration k = 1, the system now has the familiar form of equation (C.21) correspond- 
ing to the case k = 0, 


(C.47) 



i yy 

1 xy 
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A\ xx 


V i 
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B\y 
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U\ + 


L\y 

L lx. 


y i 


(C.49) 


For the above system, we apply the same procedure as developed for the case of A: = 0. The 
SCB algorithm continues until the term A k s £ is equal to zero. At the last fc £/l -iteration, the 
subsystem will be of the following form 


'Vk 


[ A kyy A kyy A kyx 1 


yk' 


't Bky 0 L k y ' 


’ Uk ' 
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In the above recursive algorithm, the system states will have the following components 


•Escb — 


y o 
v o 
V 0 
V 1 
V 1 
Vi 

Vj 
Vi 
Vi 

Vk 

Vk 

x k 

We note the following useful relations 

dim(yj) = dim{uj) 


(C.51) 


dim(yj- 1 ) = dim{ 


Vi 
Vj 
Vj J 


) 


(C.52) 


for 0 < j < k. 

Finally, we can partition x sc b further into the set of y/, y &, x c and x a states defined as 
follows, 

'ycT 
y o 
v i 

Vi 


Vf 


Vj 

Vi 

The number of infinite zeros of order j is equal to the dimension of the vector y 3 (or u 3 ) for 
0 <j<k. 

"ft)' 


(C.53) 


Vb = 


y i 

Vi 


(C.54) 


L^J 

The system structure associated with the state equations for yj has no direct input terms 
from Uj, u k and x 3 for 0 < j < k. Hence, the presence of any of these states would indicate 
that the overall system is not right-invertible. 
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If the term Bkx at the last iteration is nonzero, then the set of inputs Uk exists and the 
system will not be left-invertible. In this case, one can further partition the states Xk into 
x a and x c states where x c belongs to the controllable subspace of Uk- Namely, 


' X c 


.1 

Aca 


* x c 

+ 

0 

B c 

L c ] 


* 

_X a , 


. 0 

^4ao . 


X a _ 

.0 

0 

L a _ 


U k 

.Ilk. 


(C.55) 


Eigenvalues of the matrix A aa corresponding to the uncontrollable subsystem x a are simply 
the system invariant zeros. 
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Appendix D 


Robust Gradient Routines 

intgl2 

• FORTRAN subroutine for computing: X — [ e Fp * s CovO e Fp '* s ds 

Jo 

• INPUTS: 

Nsysord is the order of the system matrix 

N2sysord is the row dimension of CovO, usually = Nsysord 

Fp, the system matrix 

CovO, the base matrix within the integral (see above formula) 
work, an double-precision valued work array 
iwork, an integer valued work array 
tf, the final time 

• OUTPUTS: 

X matrix 

• SUBROUTINE CALLS: 
abx matrix multiplication 

abtx matrix multiplication of one matrix and the transpose of another 
atbx matrix multiplication of the transpose of a matrix and another 
abte transpose of a matrix 

lsolve solves a system of equations, often used to compute matrix inverses 


PAGE 
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conex2 


• FORTRAN subroutine for computing 

M = Jo Jo eFp '* iv ~ s)ct Q C eFp * V * Gam P * Gam P' * e Fp '* s ds dv 


• INPUTS: 

Nsysord is the order of the system matrix 

N2sysord is the row dimension of CtQC, usually equal to Nsysord 

Nnoises is the column dimension of Gamp, or the number of independent noises 

Fp, the system matrix 

CtQC, the base matrix within the integral (see above formula) 

Gamp, the secondary base matrix within the integral (see above formula) 
work, a double-precision valued work array 
iwork, an integer valued work array 
tf, the final time 

• OUTPUTS: 

DCost, the M matrix 

• SUBROUTINE CALLS: 
abx matrix multiplication 

abtx multiplication of one matrix and the transpose of another 
atbx multiplication of the transpose of one matrix and another 
abte transpose of a matrix 

lsolve solves a system of equations, often used to compute matrix inverses 
PadExp is the Pade matrix exponentiation routine (without squaring or scaling) 
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PadExp 

• FORTRAN function for matrix exponential Pad4 series specialized to this case of 
function evaluation — where the matrix has a known structure 


r 

' Fp' CtQC 0 

1 


Qn Rn Y n 


0 —Fp e Fptt Cov 0 

dT\ 

► = 

0 Sn Un 

< 

l 

O 

o 

•5 

J 


0 0 Qn 


• INPUTS: 

Nsysord is the order of the system matrix 

N2sysord is the row dimension of CtQC, usually equal to Nsysord 

serlen is the length of Pade series desired 

iwork, an integer valued work array 

dT, the final time (usually a small increment here) 

Fp, the system matrix 

CtQC, a base matrix within the M integral 

CovO, a base matrix within the M integral 

Qn,Rn,Sn,Un,Yn are all numerator summands for exponential series of a matrix 3 
times the dimension of Fp (5 submatrices, not 9, because some terms are 0 and 
others are repeated) 

Qd,Rd,Sd,Ud,Yd are all denominator subblocks of the exponential series. These 
are nominally work arrays 

workl,work2 are temporary storage work matrices (same size as Fp) 
work, a double-precision valued work array 

• OUTPUTS: 

Yn is the M. matrix for small scaling time. It is part of the matrix exponential 
Rn, a part of the overall matrix exponential ([1,2] partition above) 

Un, a part of the overall matrix exponential ([2,3] partition above) 

Sn, the exponential of the -Fp matrix and part of the overall matrix exponential 
above 

Qn, the exponential of the Fp’ matrix and part of the overall matrix exponential 
above 

• SUBROUTINE CALLS: 
abx matrix multiplication 

abtx multiplication of one matrix and the transpose of another 
atbx multiplication of the transpose of a matrix and another 
abte transpose of a matrix 

lsolve solves a system of equations, often used to compute matrix inverses 
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testintg 

• TESTINTG is a program for testing the intgl2 subroutine. This test is by way of com- 
parison with the “diagonal” method for a series of test cases, both with diagonalizable 
and non-diagonalizable matrices 

• INPUTS: 

Test inputs are a series of system matrices put in the program via FORTRAN DATA 
statements (they are basically internal to the program) 

• OUTPUTS: 

Test results that list the input matrices and the results of the computations for both 
the diagonal and robust routines are displayed directly to the screen 

• SUBROUTINE CALLS: 

testin is the actual subroutine that does the tests, but it is part of the same source 
file, so it is not all that separate 

eigenv does an eigenvalue-eigenvector decomposition of the system matrix 
abe assigns one matrix to another 

lsolve solves systems of equations, in this case it is used to find a matrix inverse 
fexp is a specialized scalar exponentiation routine particular to SANDY 
aatx multiplies a matrix by its transpose 
abx multiplies two matrices 

abtx multiplies one matrix by the transpose of another 
intglx is the original “diagonal” means of calculating X 
prntm prints a matrix out to the screen with a banner and label 
intgl2 is the “robust” subroutine for computing X 
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testcone 


• TESTCONE is a program for testing the conex2 subroutine. This test is by way of 
comparison with the “diagonal” method for a series of test cases, both with diagonal- 
izable and non diagonalizable matrices 

• INPUTS: 

Test inputs are a series of system matrices put in the program via FORTRAN DATA 
statements (they are basically internal to the program) These tests are the same 
as for the program testintg 

• OUTPUTS: 

Test results that list the input matrices and the results of the computations for both 
the diagonal and robust routines are displayed directly to the screen 

• SUBROUTINE CALLS: 

tstrun is the actual subroutine that does the tests, but it is part of the same source 
file, so it is not all that separate 

eigenv does an eigenvalue-eigenvector decomposition of the system matrix 
abe assigns one matrix to another 

lsolve solves systems of equations, in this case it is used to find a matrix inverse 
fexp is a specialized scalar exponentiation routine particular to SANDY 
aatx multiplies a matrix by its transpose 
abx multiplies two matrices 

abtx multiplies one matrix by the transpose of another 
atbx multiplies the transpose of one matrix by the other 
conexp is the original “diagonal” means of calculating M 
prntm prints a matrix out to the screen with a banner and label 
conex2 is the “robust” subroutine for computing M 
intgl2 is the “robust” subroutine for computing X 
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Appendix E 

General Utility Routines for 
Synthesis and Analysis 


pseudoi.m 

• Tinv = pseudoi(T); 

MATLAB function for computing pseudo-inverse of T via the Singular Value Decom- 
position. An inverse of T is guaranteed in the sense that Tinv*T = I’, where I’ will 
have either l’s or 0’s along the diagonal, and 0’s everywhere else. The number of l’s 
is equal to the rank of T 

• INPUTS: 

T is the input matrix 

• OUTPUTS: 

Tinv is the pseudo-inverse 


rmodal.m 

• [T,D] = rmodal(A); MATLAB routine for real eigenvalue/eigenvectors 

• INPUTS: 

A matrix to be decomposed 

• OUTPUTS: 

T real transformation matrix 

D “Diagonalized” matrix of real eigenvalues and 2x2 cr-u blocks 
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balreal2.m 


• [Abal,Bbal,Cbal,HSV,T] = balreal2(A,B,C); 

Robust form for balanced realization: algorithm works even when the system contains 
uncontrollable/unobservable modes 

• INPUTS: 

A,B,C, System dynamics, input distribution, and output distribution matrices 

• OUTPUTS: 


AbaI,Bbal,Cbal Balanced state matrices 

HSV are the Hankel singular values 

T is the transformation matrix for balanced realization 


gainloci.m 

• [zeroid] = gainloci(A,B,C,D,K, output, input, gains); 

MATLAB routine to compute closed loop eigenvalues for field of gains (negative feed- 
back only). Normally used in single loop gain margin computation 

• INPUTS: 

A,B,C,D System matrices 

K Nominal feedback gain (scale factor) for loop to be evaluated 
output is the index of the output for the loop to be fed back 
input is the index of the input of the loop fed back 
gains is the field of gains, relative to K 

• OUTPUTS: 

zeroid is the closed loop eigenvalue series corresponding to the gain field 
Gain margin in both magnitude and db 
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phasloci.m 

• [zeroid] = phasloci(A,B,C,D,K, output, input, phases); MATLAB routine to compute 
closed loop eigenvalues for a set of phase values (deg), negative feedback only. Nor- 
mally used in single-loop gain margin computation 

• INPUTS: 

A,B,C,D System matrices 

K Nominal feedback gain (scale factor) for loop to be evaluated 
output is the index of the output for the loop to be fed back 
input is the index of the input of the loop fed back 
phases is the set of phases (deg) 

• OUTPUTS: 

zeroid is a series of closed-loop eigenvalues corresponding to the given set of phases 
Phase margin in degrees 


lqrcross.m 

• [K,S] = lqrcross(A,B,C,D,Q,R); 

MATLAB routine to compute state feedback gains and Riccati solution by direct 
eigenvalue/eigenvector partitioning 

• INPUTS: 

A,B,C,D System matrices, with control input u in, criterion output z 
Q Criterion weighting matrix 
R Control weighting matrix 

• OUTPUTS: 


K State feedback matrix 
S Riccati solution 


maklticl.m 


• [Acl,Bcl,Ccl,Dcl] = maklticl(A,B,C,D,Bi,Di,Ac,Bc,Cc,Dc) 

Creates Linear Time Invariant time domain dosed loop system matrices from plant 
(A, B,C,D,Bi,Di) and controller (Ac,Bc,Cc,Dc) matrices. The plant has both a control 
input u and a command input Uc: 

dx 

— = A*x + B*u + Bi* Uc(often B=Bi, but not always) 

dt 

y = C*x + D*u + Di*Uc 


Assuming NEGATIVE feedback from the controller, this routine uses the following 
equations: 


Acl 

Bel 

Cel 

Del 


A — B(I + Dc* D)~ 1 Dc*C —B(I + Dc * D)~ l Cc 
Bc(I + D * Dc)~ l C Ac - Bc(I + D* Dc)~ l D * Cc 

Bi — B(I + Dc* D)~ l Dc * Di 

Bc(I + D * Dc )- 1 D 

(/ + D* Dc)~ l C -(I + D* Dc)~ 1 D*Cc~ 

-(I + Dc* D)~ 1 Dc*C -{I + Dc* D)- x Cc 

( I + D*Dc)- 1 Di 
— (/ + Dc* D )- 1 Dc * Di 


If either (I + D * Dc) or (I + Dc* D) are non-invertible, the solution is beyond the 
scope of this subroutine 

• INPUTS: 

A,B,C,D,Bi,Di are the plant matrices, with Bi and Di being the command input 
distribution matrices 

Ac,Bc,Cc,Dc are the controller matrices 

• OUTPUTS: 

Acl, Bel, Cel, Del are the closed loop output matrices. The columns of Bel and Del 
correspond to the command input U c . 
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sandy2.exe 


• Executable program for controller design via optimization. This form of SANDY uti- 
lizes both diagonal and robust forms of the gradient computation algorithm. Switches 
between algorithms according to the condition number of the eigenvector matrix de- 
rived from the system dynamics matrix. 

• INPUTS: 

Data file for describing the plant, controller design parameters, and current con- 
troller setup. 

• OUTPUTS: 

Log file for the optimization run, reports how well the optimization did and many 
characteristics of the final closed loop system. 

SANDY formatted file for describing current controller 

MATLAB formatted file for describing current controller 


wcsandy.exe 

• Executable program for controller design via optimization. This form of SANDY 
utilizes the diagonal form of the gradient computation algorithm. This program also 
uses the variant “worst-case” cost function approach. 

• INPUTS: 

Data file for describing the plant, controller design parameters, and current con- 
troller setup. 

• OUTPUTS: 

Log file for the optimization run, reports how well the optimization did and many 
characteristics of the final closed loop system. 

SANDY formatted file for describing current controller 
MATLAB formatted file for describing current controller 
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Appendix F 


UH— 60 Rotocraft Design Files 
(Chapter 7) 


uh60comd.m 

• MATLAB script file for generating and plotting command responses 

• INPUTS: (no formal argument list-items named actually used) 

Acl Nomenclature for dynamics matrix 
Bcmd Command input distribution matrix 

Ceval Output distribution matrix — assume the outputs are in the order: p, q, r, u, 
v, w, (j>, 6, So, Ss, Sc, and Str- 

Dcmd Command to output direct distribution 
icmd Selects appropriate column of Bcmd 
U2 Command input profile 
T2 Time base 

• OUTPUTS: 

Y, the generated responses from the MATLAB function Isim. The plots of p, q, r, 
u, v, w, (f>, 6, So, Ss, Sc, and Str are in 3 sets of panels. 

Plots of command responses. 


makideaO.m (makideal.m) 

• MATLAB script file to fabricate the ideal response model in the system matrices: 
Aideal, Bideal, Cideal, and Dideal. 

• OUTPUTS: 

Aideal, Bideal, Cideal, Dideal are the literal matrices (namewise) created-no argu- 
ment list for input or output. 11th order system. 


PAGE _ / L. INTENTIONALLY BLANK 


PMH5SWNG PAGE BLANK NOT FILMED 
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ad2delay.m 

• [Ad,Bd,Cd,Dd] = ad2delay(A,B,C,D,Tdelay); 

MATLAB function to add 4 first-order actuator delays to the system model. 

• INPUTS: 

A,B,C,D is the system model, assumed with 7 inputs. 

Tdelay is the delay time, in seconds. 

• OUTPUTS: 

Ad,Bd,Cd,Dd is the delayed actuator output system. 


ado delay, m 

• [Ad,Bd,Cd,Dd] = adodelay(A,B,C,D, Tdelay); 

MATLAB function to add 12 first order sensor delays to the system model. 

• INPUTS: 

A,B,C,D is the system model, assumed with 31 states and either 12 or more sensor 
outputs. 

Tdelay is the delay time, in seconds. 

• OUTPUTS: 

Ad,Bd,Cd,Dd is the delayed sensor output system. 
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fabsylqO.m (fabsylql.m) 

• MATLAB script file for design of state feedback for the hover (forward flight) flight 
condition. Used in Luenberger Observer design; includes delayed outputs to integra- 
tors. 

• INPUTS: 

f01_g2r2.mat (fl5_g2r2.mat) UH 60 state model linearized from nonlinear simu- 
lation program output. Flight condition point: 1 knot forward velocity (15 knots 
forward velocity) . Additional environmental parameters associated with this 
flight condition: air density 0.002377 slug/cubic foot; rotor speed 27 rad/s; and 
gross weight 16,800 lbs. 

Q and R weights are internal to the file, nonetheless they are nominal “inputs”. 

• OUTPUTS: 

State feedback matrix K in file: kmataltO.mat (kmataltl.mat). 

Plots of command responses for closed loop design 

• CALLS TO: 

makideaO.m (makideal.m) produces the ideal response/feedforward. 
ad2delay.m adds 4 actuator delays to the system model, 
adodelay.m adds 12 sensor delays to the system model. 

lqrcross.m solves for state feedback gains via partitioning the Hamiltonian matrix. 
uh60comd.m to generate the plots of the command responses 
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fabuhlqO.m (fabuhlql.m) 

• MATLAB script file for design of state feedback for the hover (forward flight) flight 
condition. Used as baseline state feedback design for comparison to output feedback 
designs. 

• INPUTS: 

f01_g2r2.mat (fl5_g2r2.mat) UH-60 state model linearized from nonlinear simu- 
lation program output. Flight condition point: 1 knot forward velocity (15 knots 
forward velocity) . Additional environmental parameters associated with this 
flight condition: air density 0.002377 slug/cubic foot; rotor speed 27 rad/s; and 
gross weight 16,800 lbs. 

Q and R weights are internal to the file, nonetheless they are nominal “inputs”. 
These are the same as those of fabsylqO (fabsylql ). 

• OUTPUTS: 

State feedback matrix K in file: kmatO.mat (kmatl.mat). 

Plots of command responses for closed loop design 

• CALLS TO: 

makideaO.m (makideal.m) produces the ideal response/feedforward. 
ad2delay.m adds 4 actuator delays to the system model, 
adodelay.m adds 12 sensor delays to the system model. 

lqrcross.m solves for state feedback gains via partitioning the Hamiltonian matrix. 
uh60comd.m to generate the plots of the command responses 


evalsysO.m (evalsysl.m) 

• MATLAB script file for generating system matrices used in the ADS-33C evaluation 
file evalhovr.m ( evalford.m ). 

• INPUTS: 

f0l_g2r2.mat (fl5_g2r2.mat) UH-60 state model linearized from nonlinear sim- 
ulation program output. 

K, the state feedback matrix from the file kmataltO.mat (kmataltl .mat). 

• OUTPUTS: 

Aeval, Ceval with outputs: p, q , r, it, v, z, (f>, 0, ip, x, y, z. 

Bgust, Dgust with inputs in the principal gust directions 

Bcmd, Dcmd with input terms in the order: pitch, roll, yaw rate, heave, yaw, 
altitude. 

G, D with standard control inputs: So, Ss, Sc, and Str. 
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fabsyloO.m (fabsylol.m) 

• This MATLAB script file does a Luenberger Observer output feedback controller 
design via CLTR. 

• INPUTS: 

P01_g2r2.mat (fl5_g2r2.mat) contains the UH-60 system model for the appro- 
priate flight condition. 

kmataltO.mat (kmataltl.mat) is the file containing the K matrix from the state 
feedback solution. 

• OUTPUTS: 

Luenberger Observer dynamics matrices: [AlCjBlCjCl^GlCjPlc]. 
fabsyloO.mat (fabsylol.mat) is the file containing these matrices. 

Plots of responses. 

• CALLS TO: 

makideaO.m (makideal.m) produces the ideal response/feedforward. 

ad2delay.m adds 4 actuator delays to the system model. 

adodelay.m adds 12 sensor delays to the system model. 

scb.m computes the Special Coordinate Basis (SCB) for the given dynamics. 

lqrcross.m solves for state feedback gains via partitioning the Hamiltonian matrix. 

uh60comd.m is a script file for plotting the responses to an input. 


fabroloO.m (fabrolol.m) 

• MATLAB script file for robust analysis of Luenberger Observer based controller gen- 
erated in fabsyloO (fabsylol). 

• INPUTS: 

fabsyloO.mat (fabsylol.mat) , the data files containing the Luenberger Observer- 
based controller. 

• OUTPUTS: 

Multiloop actuator phase/gain margins 
Single Loop actuator phase and gain margins 
Single Loop sensor phase and gain margins 

• CALLS TO: 

gainloci computes gain margins for single loop, 
phasloci computes phase margins for single loop. 
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evalosyO.m (evalosyl.m) 

• MATL AB script file to generate system dynamics matrices from Luenberger Observer- 
based controller for use in evalhovr.m ( evalford.m ). 

• INPUTS: 


fabsyloO.mat (fabsylol.mat) the source data on the Luenberger Observer based 
controller. 

• OUTPUTS: 


Aeval, Ceval with outputs: p, q, r, u, v, z, <p, 0, ip, x, y, z. 

Bgust, Dgust with inputs in the principal gust directions 

Bcmd, Dcmd with input terms in the order: pitch, roll, yaw rate, heave, yaw, 
altitude. 

G, D with standard control inputs: <5 0 , 6s, 6c, and 6tr- 


fabsyroO.m (fabsyrol.m) 

• MATLAB script file to take the balanced realization of a Luenberger Observer based 
controller and look at its command responses. Internal straps allow 10th, 14th, or 
25th order controller. 

• INPUTS: 

balrdloO.mat (balrdlol.mat) is the data file containing the balanced realization 
of the Luenberger system [Alc,Blc,Clc,Glc,Plc] 

fi01_g2r2.mat (fl5_g2r2.mat) contains the UH-60 system model for the appro- 
priate flight condition. 

• OUTPUTS: 

Plots of responses. 

• CALLS TO: 

makideaO.m (makideal.m) produces the ideal response/feedforward. 

ad2delay.m adds 4 actuator delays to the system model. 

adodelay.m adds 12 sensor delays to the system model. 

scb.m computes the Special Coordinate Basis (SCB) for the given dynamics. 

uh60comd.m is a script file for plotting the responses to an input. 
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fabroroO.m (fabrorol.m) 

• MATLAB script file to perform robust analysis of reduced order Luenberger-based 
controller. Internal straps allow 25th or 14th order. 

• INPUTS: 

balrdloO.mat (balrdlol.mat) , the data files containing the balanced Luenberger 
Observer-based controller. 

fabsyloO.mat (fabsylol.mat) , the data files containing the Luenberger Observer- 
based controller. 

• OUTPUTS: 

Multiloop actuator phase/gain margins 
Single Loop actuator phase and gain margins 
Single Loop sensor phase and gain margins 

• CALLS TO: 

gainloci computes gain margins for single loop, 
phasloci computes phase margins for single loop. 


evalrosO.m (evalrosl.m) 

• MATLAB script file to generate system dynamics matrices from a balanced and 
reduced Luenberger Observer-based controller for use in evalhovr.m ( evalford.m ). 
Straps for setting the reduction at 6th, 14th, and 25th order are available. 

• INPUTS: 


fabsyloO.mat (fabsylol.mat) the source data on the Luenberger Observer-based 
controller. 

balrdloO.mat (balrdlol.mat) the source data for a balanced Luenberger Observer- 
based controller. 

sanltreOo.m (sanltrelo.m) the source files for the 7th order reduced feedforward. 

• OUTPUTS: 

Aeval, Ceval with outputs: p, q, r, u, v , z, <p, 9, ip, x, y , z. 

Bgust, Dgust with inputs in the principal gust directions 

Bcmd, Dcmd with input terms in the order: pitch, roll, yaw rate, heave, yaw, 
altitude. 

G, D with standard control inputs: So, 6s, 6c, and 6tr- 
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sanltrfO.m (sanltrfl.m) 

• MATLAB script file to assemble reference system dynamics and controller (based on 
a 14th order observer based on a reduced Luenberger observer) set-up for SANDY 
run. This SANDY run is to minimize the white noise response error (includes shaped 
white noise) between the full-order Luenberger observer-based design and the lower 
order controller. 

• INPUTS: 

balrdloO.mat (balrdlol.mat) is the data file containing the balanced realization 
of the Luenberger system [Alc,Blc,Clc,Glc,Plc] 

sanltreOo.m (sanltrelo.m) which contain the 7th order reduced and recovered 
feedforward dynamics. 

makideaO.m (makideal.m) is the script file that produces the ideal closed-loop 
system response dynamics [Aideal,Bideal,Cideal,Dideal]. 

• OUTPUTS: 

Plots of the closed-loop response using the reference controller for verification. 

Plots of the closed-loop response using the initial guess of the reduced controller for 
verification. 

sanltrflO.mat (sanltrfl.mat) Reference and controller matrices in the plant — controller 
format acceptable to SANDY. The SANDY input file sanltrf0.dat (sanltrfl .dat) 
was eventually produced from this. 

• CALLS TO: 

rmodal.m real-valued eigenvalues and eigenvectors-complex eigenvalues show up as 
2x2 blocks. 

uh 60 comd.m is a script file for plotting the responses to an input. 
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sanltrhO.m (sanltrhl.m) 

• MATLAB script file to assemble a reference system based on state feedback on the 
UH 60 as well as open loop UH 60 dynamics. The open-loop dynamics of the refer- 
ence will be closed by an initial controller guess and the recovery error defined by the 
difference in response between the state-feedback based system and the 14th order 
observer-based controller stabilized system given the same disturbance inputs. 

• INPUTS: 

sanltrfO.m (sanltrfl.m) is the data file containing the initial controller guess (from 
the model matching run). 

sanltreOo.m (sanltrelo.m) which contain the 7th order reduced and recovered 
feedforward dynamics. 

kmataltO.mat (kmataltl.mat) is the file containing the K matrix from the state 
feedback solution. 

ft)l — g2r2.mat (fl5_g2r2.mat) contains the UH-60 system model for the appro- 
priate flight condition. 

• OUTPUTS: 

Plots of the closed-loop response using the reference controller for verification. 

Plots of the closed-loop response using the initial guess of the reduced controller for 
verification. 

sanltrhO.mat (sanltrhl.mat) Reference and controller matrices in the plant — controller 
format acceptable to SANDY. The SANDY input file sanltrh0.dat (sanltrhl.dat) 
was eventually produced from this. 

• CALLS TO: 

rmodal.m real-valued eigenvalues and eigenvectors-complex eigenvalues show up as 
2x2 blocks. 

uh60comd.m is a script file for plotting the responses to an input. 
ad2delay.m adds 4 actuator delays to the system model, 
adodelay.m adds 12 sensor delays to the system model. 


159 



fabroIrO.m (fabrolrl.m) 

• MATLAB script file to perform robust analysis of a recovered 14th order controller 
originally reduced from the Luenberger Observer design. Also strapped for possible 
use of the numerical CLTR design. 

• INPUTS: 


f01_g2r2.mat (fl5_g2r2.mat) contains the UH-60 system model for the appro- 
priate flight condition. 

sanltrfDn.m (sanltrfln.m) which contain the recovered and reduced 14th order 
controller component (sanltrhOn.m and sanltrhln.m are used given a strap set- 
ting). 

sanltreOo.m (sanltrelo.m) which contain the 7th order reduced and recovered 
feedforward dynamics. 

• OUTPUTS: 

Multiloop actuator phase/gain margins 
Single Loop actuator phase and gain margins 
Single Loop sensor phase and gain margins 

• CALLS TO: 

ad2delay.m adds 4 actuator delays to the system model, 
adodelay.m adds 12 sensor delays to the system model, 
gainloci.m computes gain margins for single loop, 
phasloci.m computes phase margins for single loop. 
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evaltrfO.m (evaltrfl.m) 

• MATLAB script file to generate system dynamics matrices from a balanced reduced 
and recovered 14th order controller based on the Luenberger Observer-based con- 
troller. Output for use in evalhovr.m (evalford.m). Does some response plotting to 
verify that a proper assembly of a closed-loop system has been made. 

• INPUTS: 

f01_g2r2.mat (fl5_g2r2.mat) contains the till 60 system model for the appro- 
priate flight condition. 

sanltrfOn.m (sanltrfln.m) which contain the recovered and reduced 14th order 
controller component. 

sanltreOo.m (snaltrelo.m) which contain the 7th order reduced and recovered 
feedforward dynamics. 

• OUTPUTS: 

Aeval, Ceval with outputs: p, q, r, u, v, z, (f>, 9, xp, x, y, z. 

Bgust, Dgust with inputs in the principal gust directions 

Bcmd, Dcmd with input terms in the order: pitch, roll, yaw rate, heave, yaw, 
altitude. 

G, D with standard control inputs: <5o, 6s, Sc, and Str- 
uh60comd.m is a script file for plotting the responses to an input. 
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evaltrhO.m (evaltrhl.m) 

• MATLAB script file to generate system dynamics matrices from a balanced reduced 
and CLTR recovered 14th order controller based on the Luenberger Observer-based 
controller. Output for use in evalhovr.m (evalford.m) . Does some response plotting 
to verify that a proper assembly of a closed-loop system has been made. 

• INPUTS: 


ft)l— g2r2.mat (fl5_g2r2.mat) contains the UH-60 system model for the appro- 
priate flight condition. 

sanltrhOn.m (sanltrhln.m) which contain the recovered and reduced 14th order 
controller component. 

sanltreOo.m (snaltrelo.m) which contain the 7th order reduced and recovered 
feedforward dynamics. 

• OUTPUTS: 

Aeval, Ceval with outputs: p, q , r, u, v, z, (p, 9, ip, x, y, z. 

Bgust, Dgust with inputs in the principal gust directions 

Bcmd, Dcmd with input terms in the order: pitch, roll, yaw rate, heave, yaw, 
altitude. 

G, D with standard control inputs: <5o, 6s, 6c, and 6tr- 

uh60comd.m is a script file for plotting the responses to an input. 
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sanl4do0.m (sanl4dol.m) 

• MATLAB script file to assemble plant and controller matrices for a direct optimization 
ran. The controller form is the same as for the recovery cases: 14th order Luenberger 
observer with integral control and 7th order feedforward dynamics. 

• INPUTS: 

sanltrfO.m (sanltrfl.m) is the data file containing the initial controller guess (from 
the model matching run). 

sanltreOo.m (sanltrelo.m) which contain the 7th order reduced and recovered 
feedforward dynamics. 

f01_g2r2.mat (fl5_g2r2.mat) contains the UH-60 system model for the appro- 
priate flight condition. 

• OUTPUTS: 

Plots of the closed-loop response using the initial guess of the controller for verifica- 
tion. 

sanl4do0.mat (sanl4dol.mat) Reference and controller matrices in the plant — 
controller format acceptable to SANDY. The SANDY input file sanl4do0.dat 
(sanl4dol.dat) was eventually produced from this. 

• CALLS TO: 

uh60comd.m is a script file for plotting the responses to an input. 
ad2delay.m adds 4 actuator delays to the system model, 
adodelay.m adds 12 sensor delays to the system model. 


T 
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fabrodoO.m (fabrodol.m) 

• MATLAB script file to perform robust analysis of a 14th order controller produced 
by direct optimization. 

• INPUTS: 

f01_g2r2.mat (fl5_ g2r2.mat) contains the UH 60 system model for the appro- 
priate flight condition. 

sanl4do0n.m (sanl4doln.m) which contain the optimized 14th order controller 
component. 

sanltreOo.m (sanltrelo.m) which contain the 7th order reduced and recovered 
feedforward dynamics. 

• OUTPUTS: 

Multiloop actuator phase/gain margins 
Single Loop actuator phase and gain margins 
Single Loop sensor phase and gain margins 

• CALLS TO: 

ad2delay.m adds 4 actuator delays to the system model, 
adodelay.m adds 12 sensor delays to the system model, 
gainloci.m computes gain margins for single loop, 
phasloci.m computes phase margins for single loop. 
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eval4do0.m (eval4dol.m) 

• MATLAB script file to generate system dynamics matrices from a direct optimized 
14th order controller based on the Luenberger Observer-based controller. Output for 
use in evalhovr.m ( evalford.m ). Does some response plotting to verify that a proper 
assembly of a closed-loop system has been made. 

• INPUTS: 

f01 — g2r2.mat (fl5_g2r2.mat) contains the UH-60 system model for the appro- 
priate flight condition. 

sanl4do0n.m (sanl4doln.m) which contain the recovered and reduced 14th order 
controller component. 

sanltreOo.m (snaltrelo.m) which contain the 7th order reduced and recovered 
feedforward dynamics. 

• OUTPUTS: 

Aeval, Ceval with outputs: p, q, r, u, v, z, <j>, 6, xp, x, y, z. 

Bgust, Dgust with inputs in the principal gust directions 

Bcmd, Dcmd with input terms in the order: pitch, roll, yaw rate, heave, yaw, 
altitude. 

G, D with standard control inputs: 6o, 6s, 6c, and 6tr- 
uh60comd.m is a script file for plotting the responses to an input. 
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evalhovr.m (evalford.m) 

• MATLAB script file for evaluating system matrices for the ADS 33C hover (forward 

flight) requirements. 

• INPUTS: 

Output from appropriate evaluation setup m file (these files start with the word 
“eval”). 

Test parameters contained within the evaluation file 

• OUTPUTS: 

Test evaluations for many of the ADS 33C requirements 

• SUPPORT ROUTINES: Because the number of routines used in these files is large 

and they are specific to them, they are listed here. Details of each will follow this list. 

bandwdth determines bandwidth from system matrices in a manner similar to, yet 
unlike classical methods. 

peakmmax evaluates transient and steady-state behavior of output trace (Y) . 

passit evaluates a point to see if it is to the right of a given curve. 

getlstor does an optimization of a given curve trace to a first order response by 
adjusting the magnitude, delay and time constant. 

getlstda returns the RMS error between a given curve trace and a trial first order 
fit. 

pulshold checks compliance to ADS-33C 3.2.6 by analyzing response to an actuator 
pulse (hover only). 

fig5_3_4 does a test based on whether a criterion (either bank angle oscillations (Fig 
5 of ADS-33C), sideslip excursion (Fig 6 of ADS-33C), or a sort of scaled sideslip 
excursion (Fig 7 of ADS-33C) is acceptable. 

fig9_3_4 checks damping ratio and stability criteria given a sigma-omega pair. 

text934 Reports in text the results from the evaluation done in fig9_3_4. 
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bandwdth.m 

• [W180,Wgain,Wphase,Tphase] = bandwdth(Acl,Bcmd,Ccmd,Dcmd,dispflag) 
MATLAB function to determine bandwidth compatible with definition in mil stan- 
dard ADS-33C, Figure 2(3.3). This function requires the closed-loop system matrix, 
the relevant command input, and the relevant controlled output (only 1 each). A 
Bode plot of the system gain and phase of X/Xcmd is made. The frequency at the 
phase=180 point is found. Wphase is the frequency at a margin of 45 deg. from that. 
Wgain is the frequency 6dB higher than the magnitude at W180. 

• INPUTS: 

Acl,Bcmd,Ccmd,Dcmd are the system matrices for bandwidth determination. The 
loop is closed, and only one input, one output. 

dispflag controls verifying the results with a plot the user can see with the 180 
frequency and the phase and gain frequencies marked. 

• OUTPUTS: 

W180 is the frequency at which the phase is 180 degrees 

Wgain is the frequency at which the magnitude is 6db up from that at W180. 

Wphase is the frequency at which the phase is 135 degrees 

Tphase is known as the phase delay, which is the slope of the phase curve between 
W180 and the frequency where the phase lags an additional 180 degrees. 


peakmmax.m 

• [Ymax,Tmax,Ymin,Tmin,Ymax2,Tmax2,Yavg] = peakmmax(Ytrace) 

Evaluates transient and steady behavior of trace. Trace assumed to be step response 
with analysis of overshoot and undershoot. Tries (presently without sophistication) 
to estimate steady state value of the trace. With all of this, the standard functions 
jmax^ and jfind£ need alternate means. With all of this, this function does not ever 
try to extrapolate to find Yavg, since that can get one into severe trouble. Tmax and 
Tmin are indices into Ytrace where the respective max and min values are located. 

• INPUTS: 

Ytrace is the step response. 

• OUTPUTS: 

Ymax is the peak maximum. 

Tmax is the index into Ytrace where the value is assumed. 

Ymin is the first minimum after Ymax. 

Tmin is the index into Ytrace where the value Ymin is assumed. 

Ymax2 is the next local maximum after Ymin. 

Tmax2 is the index into Ytrace where the value is Ymax2 assumed. 

Yavg is the estimated steady-state value. 
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passit.m 


• status = passit (curve, point) 

This function tests to see if a point is to the right of the given curve. If so status is 1, 
if not, status is 0. In curve and point , the x coordinate is first, then y. The curve is 
assumed non-closed. Linear extrapolations are done before the begin point and after 
the end points based on the segments there. 

• INPUTS: 

curve is a set of points defining a curve. This curve could be closed, but in general 
it is assumed non-closed. 

point is the x,y coordinate of a point to be tested. 

• OUTPUTS: 

status is 1 if the point is to the right of the given curve, 0 if it is not. 


pulshold.m 

• quepasa = pulshold(Ttrace,Ytrace) 

Checks compliance with 3.2.6 of ADS-33C (Handling Qualities for Military Rotocraft). 
Ytrace is a response to an actuator pulse. To pass (0 return), it is to rise, then fall 
to less than lOpeak value within 10 seconds (hence the time trace needed too). It is 
to stay below this 10% for another 30 seconds. UCE=1 needs 10% settle within only 
20 seconds, this is checked too if the original settle is not met. The return condition 
will be -1 here. It is best to supply around 60 seconds of data (only 50 absolutely 
necessary). Failure is a return condition of 1. 

• INPUTS: 

Ytrace is the response to an actuator pulse. 

Ttrace is the corresponding time trace. 

• OUTPUTS: 

quepasa is set to 0 if the trace passes the requirement, 1 if it does not. 
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getlstor.m 

• [Mag,Delay,Tconst,corrsq] = getlstor( Times, Ydata ) 

Does first order fit y(t+Delay) = Mag*(l-exp(t/Tconst)) to given data. Computes 
a sort of chi-squared sense that indicates how well the fit actually behaved. This 
procedure is presented in Figure 8(3.3) of ADS-33C. This routine calls the MATLAB 
function nelder which will in turn call getlstda, which generates data for the opti- 
mization routine. The common block: global ydata Tset Nset must be present in 
the root m-file. 

• INPUTS: 


Ydata is the function trace. 

Times is the corresponding time trace. 

• OUTPUTS: 

Mag is the magnitude of the first order function. 

Delay is the time delay to the nominal zero point for this first-order function. 
Tconst is the time constant. 

corrsq is the ratio of the squares of the fitted curve and the raw data about the 
mean. 


getlstda.m 

• error = getlstda( params ) 

Returns error from latest set of first order function parameter fit. Parameters are: 
Mag, time delay, and time constant, in that order for the params vector The common 
block: global ydata Tset Nset must be present in the root m-file. 

• INPUTS: 


params is a 3-element vector whose meanings are: Magnitude of first-order function, 
time delay, and time constant. 

• OUTPUTS: 

error is the RMS error between the data and the fitting function. 


* 
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nelder.m 


• [x, iterats] = nelder(errorfn,parmlist,tol,prnt) 

nelder is a nonlinear optimization function (the Nelder Mead Simplex algorithm) 
that tries to adjust the given parameter vector in order to minimize the output of the 
argument function errorfn. This is often a library function within MATLAB, but not 
always, hence its inclusion here. 

• INPUTS: 

errorfn is a user-defined function which uses the parmlist entry in the nelder com- 
mand line to return the associated value of the cost function. The actual argu- 
ment entry is the name of this function’s m-file (which should be the name of 
the function as well). 

parmlist is the parameter list (array). The actual meanings of the individual ele- 
ments depends on your errorfn as nelder derives the relationships of the individual 
terms to the cost. 

tol is the stopping tolerance (defaults to 1.0e-3). 

prnt will describe each step if set to a nonzero value (defaults to 0). 

• OUTPUTS: 

iterats is the number of iterations (this argument is optional), 
x is the optimized parameter list (same length as input). 


fig5_3_4.m 

• [passit] = fig5_3_4(RR2betaP,RollLag,PhOSCAV, LineSpec) 

Uses function passit to evaluate input criteria and outputs what level performance 
they represent as per Figure 5(3.4) of the ADS-33C rotocraft handling specs. Also 
evaluate as per Figure 6(3.4) and Figure 7(3.4) given different values of LineSpec (the 
raw criterion data). Will detect whether one has one curve or two (4 cols versus 2). 
For use in forward flight tests only. 

• INPUTS: 

RR2betaP is nominally the roll rate to sideslip phase angle. 

RollLag is the component of an oscillation cycle needed to reach the first maximum 
of an overshoot response in the lateral-directional oscillation. 

PhOSCAV is either the roll oscillation envelope size or the sideslip oscillation peak 
normalized by the roll oscillation peak. 

LineSpec is the specific test characteristic line using these criteria. 

• OUTPUTS: 

passit is set to 1, 2, or 3, according to what criteria are met. 
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fig9_3— 4.m 

• [result] = fig9 — 3 4( sigma, omega ) 

returns 1 if meets the all MTE specs on sigma and omega according to Figure 9(3.4) 
of ADS 33C, which examines known lateral-directional modes. Returns 2 for Level I 
(all MTE but slalom, ground attack and air combat), Level II for these. Returns 3 for 
Level III for slalom, ground attack, and air combat, Level II for all others, Returns 4 
for Level III all around, and 5 for unstable. Used in forward flight analysis only. 

• INPUTS: 

sigma, omega are the real and complex parts of the eigenvalue under consideration. 

• OUTPUTS: 

result is a simple number from 1 to 5, according to the above description. 


text934.m 

• text934( result ) 

Displays appropriate text string according to number handed it. number same series 
as in the fig9 — 3 — 4 evaluation routine. Used in forward flight analysis only. 

• INPUTS: 

result is the result from the fig9—3—4 routine. 

• OUTPUTS: 

Text to screen to describe result. 
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Appendix G 


UH— 60 Rotocraft Models 
(Chapter 7) 


f01_g2r2.mat (Hover) 

• 31 state model linearized from nonlinear simulation program output 


P, Q, r 

Body axis attitude rates (deg/s) 

U, V, w 

Body axis velocities (ft/s) 

@y) $z 

Euler angles (deg) from terrestrial to 
rotocraft axes (pitch, roll, yaw) 

x, y, z 

Inertial positions (ft) 

Po, Pd, 
Pis, Pic 

Rotor flapping rates (deg/s) 

Co, C D, 
Cl S, Cl c 

Rotor lead-lag rates (deg/s) 

Po, Pd, 
Pis, Pic, 

Rotor flapping angles (deg) 

Co, C D, 
Cl5, Cic, 

Rotor lead-lag angles (deg) 

Ao, A 15 , 
Aic 

Inflow velocities ( 1 /sec) These are normalized 
by the rotor radius (26.83ft) 


• Flight condition point: 1 knot forward velocity. Additional environmental parameters 
associated with this flight condition: air density 0.002377 slug/cubic foot; rotor speed 
27 rad/s; and gross weight 16,800 lbs 

• Dynamics matrix F 

• Control and disturbance input distribution G 

• CONTROL INPUTS (deg): So (collective), 6s (sine), Sc (cosine), and Str (tail rotor) 

• DISTURBANCE INPUTS (ft/s): V gX) V gy , and V gz , gusts from the principal terres- 
trial frame directions 


PAGE. 
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fl5_g2r2.mat (Forward Flight) 

• 31 state model linearized from nonlinear simulation program output 

• Flight condition point: 15 knot forward velocity. Additional environmental parame- 
ters associated with this flight condition: air density 0.002377 slug/cubic foot; rotor 
speed 27 rad/s; and gross weight 16,800 lbs 

• System dynamics matrix F 

• Control and disturbance input distribution G 

• All other details same as for hover model. 
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Appendix H 


JPL Large Space Structure Design 
Files (Chapter 4) 


lscllsim.m 

• MATLAB script file for running closed and open loop system dynamics through sta- 
bility test. Runs stability test mentioned in Chapter 4, the structure is pulsed for 6.4 
seconds, then allowed to free oscillate. Also gives eigenvalues of closed loop system. 

• INPUTS: 

lsclnew.mat is the data file for the open loop antenna structure. 

Iscl4m5c.m is the SANDY controller output file normally used. Its contents have 
been copied into lscllsim.m. 

• OUTPUTS: 

Eigenvalues of the closed loop system 

Plots of the rundown responses at the two hub sensors of the antenna for both open 
and closed loop. 


Iscl3m5.dat 

SANDY input file that eventually produced final controller. Contains weights, distur- 
bance profile, and starting condition details. 
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Appendix I 

JPL Large Space Structure Model 
(Chapter 4) 

lsclnew.mat 

• MATLAB binary data file for linearized and reduced antenna model 

• System matrices a, b, c, d 

• 20th order model; all modes stable, most modes lightly damped (damping 0.007 to 

0 . 01 ) 

• 6 inputs: HAl and HA10 for hub actuators; RAl, RAA, RA7 and RA10 for rib root 
actuators 

• 6 outputs: HS 1 and HS 10 for hub sensors; RSI, RS4, RS7 and RS 10 for rib root 
sensors 

• All sensors collocated with corresponding actuators 

• All actuator inputs in Newton-meters 

• All sensor outputs in radians 


PWECECHNe PAGE BLANK NOT FUMED 
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Appendix J 


One-Dimensional Rotocraft 
Design Files (Chapter 2.7.1) 


helicm0.dat, helicml.dat, helicm2.dat, helicm3.dat 

SANDY input files to run the 4 cases of initial condition versus gust response weighting. 
Weights of 0, 0.1, 1.0 and 10. on the initial conditions are used. 

helilr0.dat, helilrl.dat, helilr2.dat, helilr3.dat 

SANDY input files to run the 4 cases of loop recovery where the actuator reliance is 
lessened by injecting fictitious noise there. Weights of 0, 0.1, 1.0, and 10. are used for 
this actuator noise. Files are input files for both regular (H 2 ) and worst-case versions of 

SANDY. 


helicm.m 

• MATLAB script file for analyzing H 2 versus worst-case algorithm performance for 
a series of designs involving increasing degrees of stability augmentation over distur- 
bance rejection. 

• INPUTS: 

helild.mat is the MATLAB-compatible data file for the basic open-loop one-dimensional 
helicopter model. 

helicmOh.m, helicmlh.m, helicm2h.m, helicm3h.m are the H 2 design output 
files from their respective SANDY runs. 

helicmOw.m, helicmlw.m, helicm2w.m, helicm3w.m are the worst -case design 
output files from their respective WCSANDY runs. 

• OUTPUTS: 

RMS values of disturbance responses for displacement and pitch 
Gain margins of the closed loop system for the 8 cases. 

Phase margins of the closed loop system for the 8 cases. 

Plots of the command responses to a position (X) change. 


PMCSDtNC PAGE BLANK NOT FR WED 
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helilr.m 


• MATLAB script file for analyzing H 2 versus worst-case algorithm performance for a 
series of designs involving increasing degrees of actuator-oriented loop recovery over 
disturbance rejection. 

• INPUTS: 

helild.mat is the MATLAB-compatible data file for the basic open-loop one-dimensional 
helicopter model. 

helilrOh.m, helilrlh.m, helilr2h.m, helilr3h.m are the H 2 design output files 
from their respective SANDY runs. 

helilrOw.m, helilrlw.m, helilr2w.m, helilr3w.m are the worst-case design out- 
put files from their respective WCSANDY runs. 

• OUTPUTS: 

RMS values of disturbance responses for displacement and pitch 
Gain margins of the closed loop system for the 8 cases. 

Phase margins of the closed loop system for the 8 cases. 

Plots of the command responses to a position ( X ) change. 
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